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ABSTRACT 


The  transient  behavior  of  the  boundary  layer  on  a  flat 
plate  following  an  impulsive  acceleration  of  the  surrounding  fluid  is 
analysed.  The  effects  of  low  density  are  partially  accounted  for  by 
using  boundary  conditions  which  allow  a  slip  velocity  and  a  temperature 
jump  at  the  wall.  The  equations  governing  the  boundary  layer  are  re¬ 
written  for  a  coordinate  system  suggested  by  the  slip  boundary  condi¬ 
tions,  and  the  resulting  equations  are  solved  numerically.  Expressions 
for  the  shear  stress  and  heat  transfer  at  the  wall  are  developed  in 
terms  of  the  numerical  solution,  and  several  solutions  are  compared 
with  those  in  the  literature. 
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CHAPTER  I 
INTRODUCTION 


1.1  Statement  of  Objectives 

There  is  at  present  an  increasing  interest  in  high-altitude 
flight  where  the  air  density  is  very  low  and  the  continuum  regime  flow 
conditions  fail  to  apply.  Thus  a  theory  which  attempts  to  describe  the 
behavior  of  high-altitude  high  speed  machines  should  take  into  account 
as  many  of  the  known  characteristics  of  low-density  gases  as  possible. 
Such  characteristics  include  dissociation  and  ionization  at  high 
temperatures,  free-molecule  interactions,  and  velocity  and  temperature 
discontinuities  at  solid  boundaries.  When  evaluating  the  forces 
acting  on  an  object  subjected  to  a  transient  high-speed  flow,  viscous 
interactions  have  to  be  analysed.  Therefore  an  investigation  into  the 
characteristics  of  a  time-dependent  boundary  layer  in  a  rarefied  gas 
is  an  attractive  and  interesting  problem.  A  sharp-leading-edge  body 
will  be  considered  in  this  investigation. 

The  problem  to  be  discussed  here  is  one  in  which  the 
transient  behavior  of  the  boundary  layer  on  a  flat  plate  following  an 
impulsive  acceleration  of  the  surrounding  fluid  is  analysed.  The 
hypersonic  assumption  [1] 


=  0  , 


(i.i) 


1  ■  /<  ^  ii  jlj  1  ^  w  f  ft 
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where 


H  =  h  + 


2 ,  2 
u  +v 


is  used  throughout  the  analysis.  These  relations  are  based  on  the 
assumption  that  the  free  stream  adjusts  itself  instantaneously  to 
the  new  conditions,  whereas  the  boundary  layer  itself  is  non-stationary 
(a  similar  case  is  mentioned  by  Lagerstrom  [2]).  The  interaction  be¬ 
tween  the  viscous  and  inviscid  effects  is  considered,  and  some  of  the 
effects  of  low  density  are  accounted  for  by  using  boundary  conditions 
which  allow  a  slip  velocity  and  a  temperature  jump  at  the  wall.  Ex¬ 
pressions  are  obtained  for  the  transient  shear  stress  and  heat  transfer 
at  the  wal  1 . 

1.2  Review  of  Pertinent  Literature 

The  presence  of  a  slip  velocity  at  the  interface  between  a 
solid  and  a  gas  at  low  pressures  was  first  deduced  by  Kundt  and  War¬ 
burg  [3]  at  the  turn  of  the  century,  when  they  observed  that  the 
damping  of  an  oscillating  disk  by  the  surrounding  gas  decreased  with 
the  pressure.  Soon  afterward  Maxwell  [4]  developed  an  approximate 
analysis  of  a  monatomic  gas  adjacent  to  an  isothermal  surface,  showing 
a  relationship  between  the  slip  velocity  and  the  gradient  of  the 
tangential  flow  velocity. 


uw  ■  A*w  > 
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(1.2) 
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More  refined  calculations,  involving  approximate  solutions  to  the 
Maxwel 1 -Boltzmann  equation  in  the  vicinity  of  the  wall,  yield  the 
results  given  by  Kennard  [5], 


u  -  /u  (iy.)  + 

UW  A\v  [dy).l  4  0w  [dx\t 

w  w  w 


(1.3) 


and 


a  =  6  +  BA  (|§-) 
w  p  w  ay 


w 


1  u  (ii) 

2  Re (  w  Wlf 

w  w 


(1.4) 


The  second  term  in  equation  (1.3)  indicates  that  a  temperature  gradient 
along  a  surface  induces  a  flow  in  the  direction  of  increasing  tempera¬ 
ture.  This  phenomenon  is  known  as  "thermal  creep"  [6]  and  is  usually 
of  little  significance  when  velocity  gradients  at  the  wall  are  large. 

It  will  be  neglected  in  the  present  work. 

Equation  (1.4)  is  known  as  the  "temperature  jump"  boundary 
condition  -  a  condition  which  exists  in  rarefied  slip  flow  [5].  For 
the  normal -dens i ty  continuum  regime  the  gas  temperature  at  the  wall 
is  equal  to  the  wall  temperature.  However,  if  the  gas  is  slightly 
rarefied,  the  gas  temperature  at  the  surface  can  differ  from  the  wall 
temperature.  A  theory  describing  this  temperature  difference  was  first 
derived  by  Poisson  [7]. 

A  question  which  is  of  prime  importance  when  considering 
low-density  gases  is  the  applicability  of  the  continuum  Navier-Stokes 
equations.  Kogan  [8]  showed  that  Prandtl's  boundary  layer  formulation 


. 

- 


* 


isSSil 


is  not  only  consistent  with  boundary  conditions  (1.3)  and  (1.4),  but 
also  provides  slip  solutions  at  least  as  accurate  as  any  produced  by 
the  higher  approximations  such  as  the  Burnett  or  Thirteen  Moment 
equations.  Kogan  noted  that  "The  Burnett  equations  permit  us  to  ob¬ 
tain  a  more  accurate  description  of  the  flow  in  a  region  where  the 
Navier-Stokes  equations  are  also  applicable,  but  they  do  not  permit 
us  to  progress  into  the  region  where  the  latter  are  inapplicable". 
This  finding  agrees  with  the  conclusions  of  Schaaf  and  Chambr£  [9] 
wherein  they  state  that  the  Navier-Stokes  formulation  is  the  best 
available.  A  solution  of  the  Maxwel 1-Bol tzmann  equation,  while  being 
the  most  realistic,  is  usually  too  complex  to  be  practical. 

Several  special  problems  in  slip  flow,  such  as  the  Couette 
flow  between  flat  plates  and  concentric  cylinders,  stagnation  point 
flow,  and  the  case  of  an  incompressible  boundary  layer  [10]  have  been 
examined  during  the  last  two  decades.  In  Ref.  [11],  slip  flow  over  a 
flat  plate  is  discussed  from  another  point  of  view  by  analogy  with 
the  first  Stokes  problem  of  an  infinite  flat  plate  suddenly  set  into 
steady  motion  in  a  viscous  fluid.  This  approach  is  relatively  simple 
but  is  limited  by  questions  of  the  validity  of  the  analogy  for  slip 
flow. 

Maslen  [12]  obtained  the  solution  to  the  problem  of  slip 
flow  in  a  compressible  boundary  layer  on  a  flat  plate  by  introducing 
perturbation  of  the  known  continuum  solution  (the  Blasius  problem). 
The  boundary  layer  equations  are  expanded  in  powers  of  the  small 
parameter  e,  where 
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Me  /y" 

/■RiT 

x 


(1.5) 


to  give  two  sets  of  equations,  one  of  zero-order  in  e,  and  the  other 
containing  the  first  order  quantities.  The  same  procedure  is  followed 
with  the  slip  boundary  conditions  (1.3)  and  (1.4).  The  zero-order 
system  is  seen  to  be  the  usual  continuum  boundary-layer  system  to 
which  the  solution  is  known,  and  the  first-order  system  is  solved 
explicitly  in  terms  of  the  zero-order  quantities.  This  technique  was 
later  applied  by  Aroesty  [13]  in  order  to  introduce  slip  effects  to 
the  strong-interaction  problem  of  Li  and  Nagamatsu  [14].  Oguchi  [15] 
combined  the  linear  perturbation  and  Karman-Pohl hausen  methods  to 
obtain  a  simpler  and  more  useful  solution  to  Maslen's  problem. 

Recently  Libby  and  Chen  [16]  presented  a  unique  solution  to 
the  laminar  boundary  layer  equation  with  uniform  injection  of  mass  at 
the  plate  surface;  the  flow  in  this  case  is  non-similar.  Wazzan,  Lind, 
and  Liu  [17]  have  found  it  possible  to  apply  the  same  technique  to 
the  problem  of  slip  at  the  wall.  The  solution  is  obtained  by  a  series 
expansion  in  terms  of  a  slip  parameter  a(S),  where 

u  (S) 

f'(S,n=0)  =  -ij - =  o(S)  ;  S  =  S(x)  .  (1.6) 

e 


With  the  above  approach  a(S)  becomes  an  independent  variable  so  that 
the  results  can  be  applied  to  any  value  of  o(S)  desired.  Since  the 
method  requires  that  a(S)  be  known  a  priori,  it  is  not  suitable  for  use 
with  the  classical  slip  boundary  conditions  (1.3)  and  (1.4). 


.  :•  $ 


6 

The  method  chosen  for  the  present  work  is  a  modification  of 
the  technique  used  by  Gupta  and  Rodkiewicz  [18]  for  unsteady  non-slip 
flow  with  strong  interactions.  Because  the  flow-field  with  slip  is 
non-similar,  it  is  found  necessary  to  introduce  an  x-dependant  vari¬ 
able  5  into  the  semi-similar  (n,x)  coordinate  system  used  by  Gupta. 

The  differential  equations  thus  obtained  are  subjected  to  the  classical 
slip  boundary  conditions  and  are  solved  by  an  iterative  numerical  scheme 
due  to  Smith  and  Clutter  [19]. 
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CHAPTER  II 


EQUATIONS  GOVERNING  THE  FLOW 


2.1  Boundary  Layer  Equations 

We  consider  the  unsteady,  two-dimensional  flow  of  a  slightly 


rarefied  compressible  viscous  fluid  over  a  semi-infini te  flat  plate. 
The  coordinate  system  is  attached  to  the  plate,  with  the  origin  at  the 
leading  edge.  Using  the  usual  boundary  layer  approximations,  the  con¬ 
servation  equations  reduce  [19]  to  the  following: 

Conservation  of  Mass 


t  +  £(p“)  +^(pv)  =  0 ; 


(2.1) 


Conservation  of  Momentum 


x: 


8P  ,  8  /  9U  \  . 

8x  3y  3y'  5 


(2.2) 


(2.3) 


Conservation  of  Energy 


■ 
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These  equations  may  be  specialized  for  the  region  outside 
the  boundary  layer,  where  viscous  and  heat  conduction  effects  are 
negligible.  In  this  region  the  velocity,  temperature,  pressure  and 
density  depend  only  on  x  and  t.  Then 


_d_ 

dx 

(Peue) 

=  0  , 

(2.1a) 

dll 

dpp 

peUe 

e 

dx 

dx  * 

(2.2a) 

and 

dh 

dpp 

pe 

_ e  . 

dx 

dx 

(2.4a) 

where 

the  hypersonic 

assumption 

(i.i) 

has  been  used 

to  remove  the 

time- 

derivatives  of 

Pe  *  ue»  ^g  > 

and  Pe 

• 

» 

Thus  the  boundary  conditions 

on  equations 

(2.1)  through 

(2.4) 

are 

y  =  0: 

u(x. 

O.t)  = 

uw(x,t)  , 

v( 

x,0,t) 

=  0  , 

(2.5) 

h(x,0,t)  =  hw(x,t)  , 


and  for  an  insulated  wall, 


ah 

ay 


=  -  Pr  u. 


y=0 


3u 
w  3y 


y=0 


- 

■ 
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y  =  ye:  u(x,ye,t)  =  Ue(x)  , 


and 


h(x,ye,t)  =  hg(x)  , 


(2.6) 


2.2  Transforming  the  Boundary  Layer  Equations 

The  previous  section  has  shown  our  problem  to  be  one  which 
requires  the  solution  of  coupled  partial  differential  equations.  There 
are  five  dependent  variables  u,  v,  p,  p,  and  h,  which  are  functions  of 
the  independent  variables  x,  y,  and  t.  The  complexity  of  this  problem 
makes  a  direct  solution  of  equations  (2.1)  through  (2.4)  a  difficult 
task.  It  is  possible,  however,  to  simplify  the  problem  by  using  a 
series  of  coordinate  transformations.  The  end  result  is  a  pair  of 
coupled  partial  differential  equations  which  are  more  amenable  to  a 
numerical  solution  because  the  number  of  dependent  variables  has  been 
reduced  from  five  to  two.  Since  no  similarity  transformation  was 
found,  we  still  have  three  independent  variables. 


2. 2. a  Restricted  Dorodnitsyn-Howarth  Transformation 

The  explicit  dependence  of  the  boundary-layer  equations  on 
density  p  may  be  partially  eliminated  by  introducing  the  Dorodnitsyn- 
Howarth  transformation  as  applied  to  the  y-variable.  Thus,  when  the 
new  independent  variables 


x  = 


x  ; 


t  =  t 


(2.7) 


are  used,  the  partial  derivatives  with  respect  to  the  old  variables 
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become 


3X 


JL  +  J_ 

3x  9x  3y 


_a_  _  __p _ 3_ 

3y  "  pref  37  ’ 


3£_3 


(2.8) 


When  these  relations  are  substituted  into  equations  (2.2)  and  (2.4), 
we  obtain,  respectively. 


.[*!Uu*!U  (f£+u£  + 


3t 


3x 


3x  p 


-P_v)  %] 

ref  3y 


3P 


A  ( ,  PP  1M 
-  P 


)  , 


(2.9) 


3x  'ref  3y  Href  3y 


and 


P[(M+M|)  +  U(M+M|* 

3t  3y  1  3x  3y  dX 

=  p  JL  r  pp  111  +  pp 
pref  3y  pref^r  3y  pref^r 


+  v 
pref  3y 

(Pr-1)  u  . 

9y 


(2.10) 


3P 

Equations  (1.1)  and  (2.3)  together  imply  that  ^  =  0  . 

The  derivation  of  equation  (2.10)  requires  the  use  of 
equations  (2.2)  and  (2.3)  in  (2.4).  The  relation 


, 
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h  =  H  -  \  , 

which  is  a  consequence  of  the  hypersonic  assumption,  is  also  used  in 
the  derivation. 

At  this  point  a  stream  function  ^  is  introduced  such  that 


From  equation  (2.7)  we  have 


P  = 


P 


& t 
ref  ay 


From  (2.8)  and  (2.12) , 


pu 


dip 

pref  3y  ' 


(2.12) 


(2.13) 


(2.14) 


This  relation  plus  equations  (2.8)  and  (2.13)  are  substituted 
into  the  continuity  equation  (2.1).  Then  (2.1)  may  be  integrated  to 
yield 


v 


ref  /ii  + 

p  ax  ay  3x  9t 


(2.15) 


if  the  assumption  is  made  that  both  y  and  ^  are  continuous  functions 
which  are  continuously  differentiable  at  least  to  second  order  [20]. 
Boundary  condition  (2.5)  is  used  to  show  that  the  constant  of  integration 
is  zero. 


If  ip  is  defined  by  equations  (2.12)  and  (2.15),  then  the 
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continuity  equation  (2.1)  is  satisfied.  Use  of  the  stream  function 
causes  the  momentum  and  energy  equations  to  become,  respectively, 


2  2  2 

3  ip  +  B^  3  ip  _ 

ByBtf  By  ByBx"  Bx  By2 

_  1  9P_+  _1 _ B_ 

p  Bx  pref  By 


pref  3y 


(2.16) 


and 

9H  +  9H  _  9H_  9^ 

9t  By  Bx  By  Bx 

*  ^4  [7^^=+  (Pr-1)  ^4  •  (2-17) 

pref  9y  pref  r  9y  pref  r  By  By 

2.2b  Howarth-Stewartson  Transformation 

The  explicit  appearance  of  p  in  equations  (2.16)  and  (2.17) 
may  be  removed  by  further  transformations.  The  new  independent  vari¬ 
ables  are 


C(4-)  (p-M  dl 

aref  Kref 


(2.18) 


Derivatives  can  be  transformed  via 


-r=  c(4~)  (4~)  -3+  . 

9x  ref  ref  9X  9x  9Y 


_3_  =  ae  J_ 
By  aref  BY  * 


(2.19) 


. 


* 


where  C,  which  has  the  form  of  the  Chapman-Rubesin  constant  [21],  is 
in  this  case  not  a  constant,  but  is  a  function  of  x  appearing  in  the 
linear  viscosity  law, 


C  e  C(x)  = 


ye 


ref 


y  *e  ‘ 

Mref 


(2.20) 


These  relations,  when  used  in  equations  (2.16)  and  (2.17), 


yield 


c(-V(4~>  ^4:+  c(4-)3(4->  -3-4(4 


ref  3Y3T 


aref'  'Pref  3Y  3X3Y 


aref'  'Pref'  ae  3X  3Y 


c(A-)3(4_)^4=.£(A.)(4_)^ 

aref  ref  3X  3Y*  p  aref  rref  3X 


+  C(A-)(b^-) 


y 


_  ref  3°^ 

aref”Pref/  pref  3Y3 


(2.21) 


and 


3H_ 

3T 


2  P  . .  a^  2  P 

+  C(-5-)  (p-M  -  C(-5-)  (p-^-) 

aref  ref  3Y  3X  aref  ref 


3^  3H_ 
3X"  3Y 


1 

pref 


(4->2  -3  + 

aref  3Y  pref  r  3Y 


yp 


(Pr-1)( 


.  (2.22) 

3Y  3 r 


The  equation  of  state  for  a  perfect  gas. 
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P  =  pR6  , 


(2.23) 


was  required  in  the  derivation  of  (2.21). 


2.2c  Temporal  Transformation 

Gupta  and  Rodkiewicz  [18]  found  that  the  conservation 
equations  (2.21)  and  (2.22)  could  be  simplified  by  a  distortion  of 
the  time  variable.  If  the  transformation 


T  = 


aQ  2  PQ 

Cf-5-)  (d-5-)  dT  , 


0 


ref 


ref 


(2.24) 


is  introduced,  then 


_3_  _  _3_  .  91  JL 

3X  3X  3X  3T 


_3_  _  _3_ 
37  '  3Y  ’ 


(2.25) 


Thus  equations  (2.21)  and  (2.22)  become,  respectively. 


32i p  dip  d2ip  dip  d2\p  H  2  M  dMe 

3Y3T  3Y  3Y3X  "  3X  y2  "  He  ref  1  e  dX 

3T  fdip  d2\p  dip  d2ip\  _  d3ip 

3*3Y  9Y9T'3T  9Y2j_Vf  3Y3 


+ 


(2.26) 
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and 


3T  3Y  3X  '  3X  3Y  „v  l3Y  3T  "  3T  3Yj 


vref  3  r3H 
Pr  3Y  L3Y 


(Pr-1)( 


3£  3i)h 

3Y  3Y2J 


(2.27) 


where  equation  (2.11)  plus  the  following  relations  have  been  used: 


ae2  =  (y-D  he  . 


* 


e  e 


jL^e 

he  3X 


5 


(y-1)  Me2  , 


1  3he 
he  3X 


2 

aref 


dM 


M 


e  dX 


(2.28) 

(2.29) 

(2.30) 

(2.31) 

(2.32) 


Equations  (2.26)  and  (2.27)  were  obtained  by  Gupta  and 
Rodkiewicz  [18]  and  are  the  flow  governing  equations  in  terms  of  the 
independent  variables  X,  Y,  and  T.  These  new  variables  are  related 
to  the  physical  variables  x,  y,  and  t  by  the  following  expressions: 


10  ;  *  WB 


■ 
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x 


y 


x  = 


c( 


aref 


o 


(2.33) 


2. 2d  Final  Transformation 

Gupta  was  able  to  further  simplify  equations  (2.26)  and 
(2.27)  by  applying  a  similarity  transformation.  In  the  present  problem, 
because  of  the  slip  boundary  conditions  at  the  wall,  the  velocity  and 
temperature  of  the  gas  adjacent  to  the  wall  do  not  remain  constant, 
but  instead  are  functions  of  x  and  t.  In  this  case  a  similarity  trans¬ 
formation  of  the  form  found  by  Rodkiewicz  and  Reshotko  [1]  and  used  by 
Gupta  [22]  leads  to  boundary  conditions  which  are  functionally  in¬ 
consistent;  the  x-  and  t-dependence  can  not  be  eliminated.  In  Appendix 
A  the  boundary  conditions  are  examined,  with  the  result  that  the 
following  coordinate  transformations  are  suggested: 


where 


L  =  L(X,T)  = 


(2.35) 


Vefprefaref 


and 


(2.36) 


■ 
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We  let 


S  =  S(£sti,t)  =  p- 

ne 


(2.37) 


The  derivatives  appearing  in  (2.26)  and  (2.27)  are  then 
transformed  as  follows: 


JL  =  MJL+iaJL+JliL 

9X  9X  9^  9X  9ri  9X  9t  * 


JL  -  in  JL 

9Y  9Y  9n  * 


_9_  =  _3_  9r  J_ 

9T  9T  9£  9T  9t  * 


(2.38) 


We  now  define  a  new  stream  function  f  such  that 


*  =  [' 


2WV/2 
7 i+ry 


•]  f(£»n,T)  . 


(2.39) 


When  equations  (2.34)  and  (2.39)  are  used  in  (2.38),  the  following  ex¬ 
pressions  result: 


LiL  =  V  f  • 
9Y  ve  T 


(2.40a) 


f ' 


(2.40b) 


' 
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9X 


Tiny 


\  Yf’ 

e 

nr ' 

■  J  Yf' 

1/2 

1/2 

X 

Ve 

f(-i  + 

_e_ 

X 


cW 
t  e 


4 


*2vref  Ve 


1/2 


Tm+ry  xi/2 


f  +  i 
T  2 


2vref  X1/2 


Tirry  v  1/2 


dV( 

dX 


m+-l 


2^7  V  1X2 
ref  e 


+  “2"  yWTf  7172“  ’ 


d2ip  _  f .  dVe  1  / (m+1 )'  Ve 

9X9Y  "  T  dX  2  v  2vpef  xl/2 


1/2 


Yf 


dve 
I  I  _ c_ 

dX 


3/2 


521  \  Kt*' 


9 


=  /5+ji 

,w2  v  2v.  ^ 


9Y 


ref  X 


3/2 

VT 


I  I 


!  1  _  (m+1 )  , 

3Y3  ~  2vref  X 


V  3/2 

FT 


9 


9H  _  „ 

9T  he  5  X  * 


(2.40c) 

(2.40d) 

(2.40e) 

(2.40f) 

(2.40g) 

(2.40h) 

(2.40i) 


. 


* 


' 


( 2 . 40 j ) 


3H  _ 
3X  " 


/istH 


2v 


ref 


u  c./l  _V _ 

e  (2  dX 


dVp  1  v1/2 

e  1  Y  e 

"  2  Y  v3/2 


) 


•  T  t  9Vo  3H 

^e^I+V^*5  af 


+  H. 


m+1 


L.  c* 
X  5 


/ - '  V  1/2 

9H  _  / (m+1 )  ve  ,, 

3Y  '  2Vf  X1'2  e 


2  V 

8  H  m+1  j.  e  r , , 

3Y2  -  2vref  He  X 


(2.40k) 


(2.401) 


If  equations  (2.34)  through  (2.40)  are  substituted  into  the 
momentum  and  energy  equations,  then  the  result  is 


2vref  ,  2vref  X  dVe  f2, 

m+1  T  m+1  V  dX  T 

e 


2vref  X  dVe  H 

m+1  V  dX  H 
e  e 


2vref  ,1.1  X  dVe.  ,, 

m+1  l2  2  V  dX  j  TT 
e 


• .  2vref 


l  i  dVo 

_  V  Tf- _ 2 _ — 

m+1  VLX  V.  dX 


+  ^  R7“  -Vff" 


dM 


(2.41) 


' 
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and 


2vref  s  +  2Vf  -  dv- 


m+1 


m+1 


T  m  _  _ §.  +  ( Ry.~ 1  \ 

T  Ll  V  dX  vi  ' 


a  dM  .  . 

rNed fms'f-sr) 


2vref  ,1  .  1  X  dVeK,f  vref  Pr-1  ae  ,,  2,-2, ,  vref  r , , 

'  “iffTT  *1  Y~  B)T)S  f - 9-e^u—  (f  +f  f  )- -ft— S 

e 


a  2  Pr  H  e 
aref  e 


Pr 


2v 


5(f*S‘-f'S*) 


(2.42) 


Now  let 


V  =  eX 
e 


m 


(2.43) 


and 


3  = 


2m 

m+1 


(2.44) 


Then,  since  H  is  independent  of  X,  and  since 

C 


Ue2  (Y-l)  Me2 
H~  “ 


1  +  M  2  ' 

2  e 


(2.45) 


equations  (2.41)  and  (2.42)  may  be  written  as 


(2-B)f‘  +  S(f '  -S)  -  2(l-gh(f'f'-ff") 


2v  1  (Y-l)  M  2  .  . 

e(fV-)  T{ - v  1  2}  (f’f'-ff")  -  ff"  -  f 

Y"‘  1  +  V-  M 

2  e 


=  5(f*f'  '-f'f*') 


(2.46) 
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and 


S  Pr[(2-3)-2(l-3)Tf']  -  S"  -  S'Pr[f-2(l-3)xf] 


9  i  (y-1 ) 

Prr  {- 


_■!  c  2)  (f'S-fS')-(Pr-l)  {- 


(y-1)  M, 


1  +  ^  M/'  ’  ’  '  1  +  ill  |*2 

2  e  2  e 


x  (f  ,2+f'f'  ")  =  Pr  £(f*S'-f'S*)  . 


■} 


(2.47) 


Since  the  use  of  the  hypersonic  assumption  limits  us  to  the 
consideration  of  large  Mach  number  flows  only,  we  may  safely  use  the 
approximation 


(y-1)  m/ 

1  +1r"ez 


}  *  2  . 


(2.48) 


This  allows  us  to  write  equations  (2.46)  and  (2.47)  as 


•  • 


(2-3) f'  +  3( i  -S)  -  2T[(^T)3+l](f,f,-ff") 


-  ff"  -  f'"  =  £(f*f"-f'f*') 


(2.49) 


and 


(2-3)  Pr  S  -  2Pr  x[(^y)3+l ](f  'S-fS ' )  -  PrfS' 


-  S"  -  2  (Pr-1 )  (f '  ,2+f  ‘  f 1  1 ' )  =  Pr  £(f*S'-f'S*)  .  (2.50) 


These  two  expressions  are  the  final  form  of  the  governing  equations 


* 


■ 


■ 
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for  the  unsteady,  two-dimensional  high-speed  flow  of  a  slightly  rare¬ 
fied  compressible  viscous  fluid  over  a  semi-infinite  flat  plate. 

The  energy  equation,  when  specialized  for  Pr  =  1 ,  becomes 

( 2-3) S  -  2T[(^y)3+l](f'S-fS')  -  fS '  -  S" 


=  S(f*S'-f'S*)  .  (2.51) 

2.3  Specialization  for  Zero  Pressure  Gradient 

Equations  (2.49)  and  (2.51)  are  valid  for  any  pressure 
gradient,  but  since  the  value  of  3  is  assumed  to  be  independent  of  x, 
the  analysis  must  be  limited  to  flows  in  which  the  pressure  gradient 
parameter  3  is  a  constant.  There  are  only  two  types  of  flow  which  meet 
this  requirement  -  the  strong  interaction  case,  for  which  we  must  have 
3  =(y-l)/y  [23],  and  the  weak  interaction  case,  with  3  =  0  .  It  is  known 
that  close  to  the  leading  edge  of  the  plate,  the  interaction  between 
the  viscous  and  inviscid  flow  regions  has  a  strong  effect  on  the  flow 
properties.  The  interaction  effects  become  weaker  as  the  distance  from 
the  leading  edge  increases.  We  thus  see  that  the  two  values  of  3  given 
above  represent  only  the  extremes  of  viscous  interaction.  There  must 
exist  a  region  on  the  plate  for  which  the  proper  value  of  3  lies  be¬ 
tween  0  and(y-l)/y  .  Since  the  formulation  of  the  problem  does  not 
allow  us  to  give  3  any  sort  of  dependence  on  x,  we  are  forced  to  deal 
with  one  of  the  two  extremes.  The  weak  interaction  case  (3=0)  is  used 
in  the  remainder  of  this  work.  There  are  two  reasons  for  this  choice: 


'  ■ 


■ 

' 
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1.  Strong  interactions  are  important  in  a  region  of  only 
vaguely  defined  size.  Although  it  is  possible  for  the 
flow  properties  and  wing  geometry  to  be  such  that  strong 
interactions  exist  over  the  entire  wing  surface,  it  is  more 
likely  that  a  wing  producing  useful  lift  would  have  a  chord 
long  enough  to  allow  an  appreciable  weakening  of  the  viscous 
interactions.  At  speeds  in  the  low  hypersonic  range  (Mach  5 
to  Mach  15)  the  region  of  strong  interactions  should  be  but 
a  small  percentage  of  the  total  wing  length.  The  observa¬ 
tions  of  the  authors  of  Ref.  [24]  indicate  that  the  strong 
interaction  region  should  be  important  only  for  flight 

Mach  numbers  greater  than  about  15. 

2.  Solving  the  coupled  strong  interaction  equations  requires 
computations  which  are  almost  seven  times  as  lengthy  as  for 
the  weak  interaction  case. 

Equations  (2.49)  and  (2.51),  when  specialized  for  zero 
pressure  gradient,  have  the  form 

2f 1  -  2T(f,f,-ff,,)  -  ff"  -  f"'  =  ^(f*f"-f'f*')  (2.52) 


and 


2S  -  2x(f 'S-fS1 )  -  fS '  -  S"  =  £(f*S  -f'S*)  .  (2.53) 


The  solution  of  these  two  equations  is  described  in  the  following  chapters. 


, 


, 


' 
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We  should  take  note  of  the  fact  that  these  equations  have 


been  developed  using  an  arbitrary  reference  condition.  The  only  as¬ 
sumption  which  has  been  made  is  that  the  reference  values  are  in¬ 
dependent  of  x.  If  no  longitudinal  pressure  gradient  exists  (3=0) 
the  free  stream  Mach  number  is  independent  of  x,  and  thus  the  free 
stream  properties  are  also  constant  and  may  be  used  as  the  reference 
values.  In  the  presence  of  a  longitudinal  pressure  gradient,  the 
free  stream  Mach  number  is  a  function  of  x,  as  are  the  free  stream 
properties.  However,  the  free  stream  stagnation  properties  remain 
constant,  and  hence  are  ideal  for  use  as  the  reference  values. 

2.4  Boundary  and  Initial  Conditions 

The  boundary  conditions  appropriate  to  equations  (2.52)  and 
(2.53)  are  as  follows: 

1.  When  y  =  ye  we  have  Y  =  Yg,  n  =  ne,  and  u  =  Ue3-  Then 


e3 


(2.54) 


Also,  we  have  H  =  H  ,  from  which 


S(£,ne»T)  =  1  . 


(2.55) 


^  M 

2.  When  y  =  0,  we  have  Y  =  0,  n  =  0,  and  u  =  AAw  .  Thus, 
from  Appendix  A, 


f'(C,0,T)  =  |  f"(?,0,T)  , 


(2.56) 


. 


. 
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so  long  as  £  >  $£  and  x  f  0.  Here  5E,  is  that  value  of  £ 
for  which  f'(S£,0,x)  =  ^- f " (6£,0,x)  =  1. 

These  restrictions  are  necessary  because  of  the  unrealistic 
behavior  of  the  boundary  condition  at  £  =  0  and  the  incon¬ 
sistency  with  condition  (2.60)  at  t  =  0. 

When  y  =  0,  we  have  Y  =  0,  n  =  0,  and  ^  =  0.  Hence 


f(S.O.i)  =  0. 


(2.57) 


At  the  plate  surface  (y=0)  we  have  6  =  0p  +  BAW 


9y 


Then,  from  Appendix  A, 


w 


S(5,0,t)  =  +  |  S'(5,0,t)  +  f"2(c,0,x)  (2.58) 

0O  5  C 


so  long  as  £  >  S£  and  x  f  0. 

For  the  case  of  no  wall  heat  transfer.  Appendix  A  shows 
that 


S'(?,0,t)  =  0. 
ad 


(2.59) 


3.  When  t  =  0,  we  have  x  =  0,  u  =  u-j  +  (Ue3-Ue-|)*,  and  hence 

U  ,  U  q 

f'(C»n»0)  =  fiUi.ni)  ri— +  (1  -  it1)  . 

1  1  1  ue3  ue3 


★ 

U-j  =  u(C,n,T<0);  u2  =  u(C,n,T=0);  u3  =  u(£,n,x>0) 


(2.60) 


. 


■ 


. 
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Also  at  t  =  0,  we  have  H  =  H^.  Then 

H2 

SU.n.O)  =  r±-  .  (2.61) 

he3 

4.  The  assumption  is  made  that  free-stream  conditions  exist 
at  the  leading  edge  of  the  plate.  Thus,  at  C  = 

f'(6£,n,x)  =  1,  S(6£,n,x)  =  1  .  (2.62) 

2.5  Correlation  with  Equations  in  the  Literature 

When  equations  (2.52)  and  (2.53)  are  specialized  for  steady 
state,  the  result  is 


f"  +  ff,.  +  £(f*f"-f'f*')  =  o  (2.63) 


and 


S"  +  fS *  +  £(f*S'-f'S*)  =  0  .  (2.64) 

Equation  (2.63)  has  been  obtained  by  Wazzan,  Lind,  and  Liu  [17] 
in  a  slightly  different  form.  A  coupled  form  of  equations  (2.63)  and 
(2.64)  was  obtained  and  solved  by  G.  Berard  [25].  A  comparison  with 
his  experimental  results  is  included  in  Chapter  V. 

For  the  case  where  there  is  no  slip  or  temperature  jump  at 
the  wall,  f  and  S  do  not  depend  on  £.  Thus  the  derivatives  denoted  by 


■ 
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are  non-existent,  and  equations  (2.52)  and  (2.53)  become 

2f '  -  2T(f'f,-ff")  -  ff"  -  f"  '  =  0 

and 

2S  -  2x(f 'S-fS' )  -  f S '  -  S "  =  0  . 

If  we  now  define  a  new  function 

g(n,x)  =  (1  +  M02)  S(n,x)  -  ^  Me2f,2(n,x) 

and  substitute  it  into  (2.66),  we  find  that 

2g(  1  -xf 1 )  -  g"  -  g'(f-2-rf)  -  (y-1)  Me2f"2 

+  (y-1)  Me2f'{2f'  +  2T(ff"-ff')  -  ff"  -  f'"}  = 

Using  (2.65)  in  this  expression,  we  have 

2g(l-Tf')  -  g"  -  g'(f-2if)-(Y-l)  Me2f  '2  =  0  . 

Equations  (2.65)  and  (2.69)  have  been  obtained  and  solved  by 
and  Reshotko  [1]. 


(2.65) 


(2.66) 


(2.67) 


(2.68) 


(2.69) 

Rodkiewicz 


Since  g,  the  non-dimensional  local  enthalpy  ratio,  is  given  by 
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9  =  i=  (1  +1^Me2)>  (2-70) 

we  may  use  (2.67)  to  obtain  an  expression  for  static  temperature  0 
referred  to  the  free  stream  stagnation  temperature  0q, 

4-  =  S(c,n,x)  -  f,2U,n,T)  .  (2.71) 

eo 

Relation  (2.71)  is  valid  only  if  the  Mach  number  is  large. 
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CHAPTER  III 

SOLUTION  OF  THE  GOVERNING  EQUATIONS 

3.1  Finding  a  Solution 

The  problem  formulated  in  Chapter  II  requires  the  solu¬ 
tion  of  a  highly  nonlinear  third  order  parabolic  system  of  partial 
differential  equations  in  which  there  are  three  independent  variables. 
Such  a  complex  problem  does  not  willingly  yield  a  closed  form  solu¬ 
tion,  so  a  suitable  numerical  technique  must  be  looked  for. 

3.2  The  Numerical  Procedure 

The  method  used  to  solve  equations  (2.52)  and  (2.53)  is 
based  on  ideas  originated  by  Hartree  and  Womersley  [26]  and  developed 
by  Smith  and  Clutter  [27]  for  the  solution  of  the  incompressible 
laminar  boundary  layer  equations.  Later  investigations  extended  the 
technique  to  the  compressible  laminar  boundary  layer  [28].  The 
method  was  also  found  capable  of  handling  binary  dissociating  mixtures 
and  transverse  curvature  effects[29  ,  30,31 ].  Many  of  the  details  of 
the  procedure  can  be  found  in  these  references. 

3.2a  Description  of  Method 

The  fundamental  idea  for  the  method  of  solution  is  that  of 
replacing  the  streamwise  ^-derivatives,  which  are  only  of  first  order, 
by  forward  finite  difference  approximations.  The  unsteady  x-deri vatives , 


. 
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also  of  first  order,  are  handled  the  same  way.  The  remainder  of  the 
equation  is  left  unchanged  and  as  a  consequence  the  equation  is  con¬ 
verted  to  an  ordinary  differential  equation  which  is  solved  over  and 
over  again  as  one  proceeds  along  the  and  x-di recti ons .  The  use  of 
forward  differences  requires  that  solutions  be  known  at  upstream 
locations. 

While  marching  in  the  ^-direction,  the  value  of  x  is  held 
constant.  When  £  has  become  sufficiently  large,  x  is  increased  by  Ax, 
the  unsteady  derivatives  are  correspondingly  updated,  and  the  marching 
procedure  is  repeated  starting  at  E,  =  0.  This  sequence  is  continued 
until  solutions  have  been  obtained  over  the  desired  range  of  £  and  x, 
e.g.  for  0  £  £  £  1 0  and  0  £  x  £  1 . 

3.2b  Modifying  the  Equations  for  Numerical  Solution 

During  the  development  of  the  numerical  method.  Smith  and 
Clutter  found  that  roundoff  errors  arising  during  the  integrations 
could  be  reduced  by  a  redefinition  of  the  dependant  variables.  We 
introduce 


F  =  f  -  n  , 

F*  =  f  -  1  , 

F"  =  f"  , 

F."  =  f...  j 


(3.1) 


■ 
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and 


G  =  S  -  1  , 

G'  =  S'  ,  (3.2) 

G"  =  S"  . 


These  variables  go  to  zero  at  the  edge  of  the  boundary  layer  whereas 
the  original  dependant  variables  f'  and  S  go  to  unity. 

By  using  the  Euler  transformation  [21]  as  applied  to  the  time 

variable. 


C  =  ,  (3.3) 

we  map  the  region  0  £  t  <_  °°  into  the  more  easily  handled  region 
0  <  c  <  1;  the  result  is  an  improvement  in  the  convergence  rate  of 
the  numerical  method. 

Applying  expressions  (3.1),  (3.2),  and  (3.3)  to  equations 
(2.52)  and  (2.53),  we  obtain  equations  which  are  well  suited  to  a 
numerical  solution.  After  slight  rearranging,  these  equations  become 


F' "  =  -  ( F+n ) F 1 '  -  C[F"  ft-  (F'+l )  ftl] 


+  2<1-?)2  If1-  2c(1-c)[(F'+1)  F"  ft] 


3? 


(3.4) 


' 


' 
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and 


G"  =  -  ( F+n ) G ’  -  C[G'  ||  -  ( F'+l )  ||] 

+  2(l-c)2  §§  -  2C(l-c)t(F'+l)  §|  -  G'  §|]  .  (3.5) 

3.2c  The  Associated  Boundary  Conditions 

It  is  of  course  necessary  to  modify  the  boundary  and  initial 
conditions  in  accordance  with  the  previous  section.  Hence  the  re¬ 
lations  (3.1),  (3.2),  and  (3.3)  are  used  in  expressions  (2.54)  through 
(2.62)  to  yield 

1 .  At  n  =  ne  , 


F'(5,ne.c)  =  0  , 

(3.6) 

G(5,ne,d  =  0  . 

(3.7) 

2. 


At  r)  =  0  , 

F'(5,0,C)  =  |  F"(?,0,?)-1  , 

F(C.O.C)  -  0 

9p  D 

G(C,0,c)  =  l  G' 

60  5 

+  l^MBLF112Uf0tC).1  > 

5 


5  >  65  .  C  1  0  ; 

(5,0,  C) 

5  >  65  ,  C  t  o  • 


(3.8) 

(3.9) 


(3.10) 


* 
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The  restrictions  on  (3.8)  and  (3.10)  are  required  in  order  to  eliminate 
inconsistencies  in  the  boundary  conditions  at  £  =  0  and  c  =  0. 

For  the  case  of  no  wall  heat  transfer,  we  have 


G‘  (£,0,c)acj  =  0  . 

(3.11) 

3.  At  C  =  0  , 

u  i 

F'  (5,n»0)  =  F-j '  (51  ,n-,)  y 

03 

(3.12) 

and 

H2 

G(5.n.O)  =  jr-  -  1  . 

He3 

(3.13) 

We  further  specify  [32]  H2  such  that 

H  (5.1,0)  =  0  . 

(3.14) 

4.  At  £  =  0  , 

F'(0,n,T)  =  0  , 

(3.15) 

G(0,ti,t)  =  0  . 

(3.16) 

Thus  for  the  third  order  momentum  equation  we  have  three 
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conditions  on  n  and  one  each  on  £  and  c.  The  energy  equation,  which 
is  of  second  order,  has  two  conditions  on  o»  one  on  £,  and  one  on  c* 

The  problem  is  therefore  completely  specified. 

3.3  The  Finite  Difference  Formulation 

The  general  form  of  the  forward  finite  difference  approximation 
to  the  £-  and  c-derivati  ves  is  as  follows: 

When  only  one  upstream  solution  is  known,  derivatives  are 
replaced  by  the  two-point  formula 


3D 
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1  =  1 


D-D  -I 
n  n- 1 

I  -I 
n  n- 1 


(3.17) 


n 


If  two  or  more  upstream  solutions  are  known,  a  more  accurate 
three-point  formula  is  used. 


3D 
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1  =  1. 


CI  -!  ,  +  I  -I  Dn 

n  n- I  n  n- ^ 


n 


'  [<In-In-l>"VrIn-2>]  ^ 


.  r  _ n  1 _  1  n 

"-2  ’ 


(3.18) 


These  are  general  equations  in  which  Dn  is  the  value  of  the  dependent 

■hh 

variable  D  at  the  ntn  I -location  In,  where  I  is  the  independant  vari¬ 
able  of  differentiation.  The  above  forms  do  not  require  a  constant 


. 
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stepsize  ini. 

3.3a  Errors  in  the  Finite-Difference  Approximations 

The  errors  introduced  by  using  expressions  (3.17)  and  (3.18) 
are  [27],  respectively 


E2  ■ 


I  -I  t 
n  n- 1 


2 

a^D, 

31^1 


(3.19) 


1  =  1 


n 


and 


E3  = 


^ In“In-l ^In"In-2^  93D 


ar 


(3.20) 


1=1 


n 


With  the  Hartree-Womersley  technique  the  most  important  error  parameter 
is  not  the  stepsize  AI  ,  but  is  instead  the  ratio  I  /AI  .  Smith  and 
Clutter  have  devoted  much  effort  to  studying  the  effects  of  this  para¬ 
meter,  and  they  have  learned  that  sensitivity  problems  and  error  pro¬ 
pagation  is  reduced  if  the  value  of  I n/AI n  does  not  exceed  25.  This 
means  that  to  have  the  same  accuracy  in  the  solution  at  all  stations  the 
stepsize  Aln  must  increase  with  increasing  In.  The  ability  to  change 
the  stepsize  AI  in  a  prescribed  manner  has  been  built  into  the  computer 


program. 


3.3b  The  Approximated  Momentum  Equation 

There  are  six  forms  of  the  momentum  equation  which  are  needed 
for  the  numerical  solution.  Each  form  is  used  at  a  different  point  in 
the  (£,c)  plane.  Figure  1  illustrates  their  use. 


. 
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Used  In  The  Numerical  Solution 

If  we  assume  that  Aln  =  ( I n-I n_i )  =  then  ^e  equations  may 

be  written  as  follows: 

Equation  ( 2) F:  This  is  a  steady-state  equation  with  2-point  derivatives 
in  and  is  used  only  at  the  points  (£i=1>  n,  Ck<Q).  Tt  has  the  form 


Fi;k  ■  -  <Fi,k+^i;k  -  e1tq:k(F1ik-F1.1,k)-(q,k+i)(qik-F1'.lik)]/A5l. 

(3.21) 


Equation  (3) F:  This  is  also  a  steady-state  equation,  but  it  uses  3- 
point  derivatives  in  £.  Used  only  at  the  points  (£.>2>n.Ck<0) ,  it  has 


the  form 


. 


1 
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Fi\k  =  "  (Fi,k+n)Fi\k  '  Ci[Filk(!  Fi,k“2Fi-l,k+  \  Fi-2sk) 


(pi,k+1^2  Fi  ,k~2Fi-l  ,k+  2  Fi-2,k^/A^i  * 


(3.22) 


Equation  ( 2 , 2 ) F :  This  is  the  first  of  the  unsteady  equations.  It 
employs  2-point  derivatives  in  both  £  and  and  is  used  only  at  the 
points  (£.=1  »n,Ck=1)  . 

F!'^  -  R.H.S.  of  (2)F.  =  2(l-?k)2(Fi>k-F->k.1)/A?k 


(3.23) 


-  qyq.k-q.k-i)^  • 

Equation  (3,2)F:  Here  we  use  a  3-point  derivative  in  £  and  a  2-point 
derivative  in  c  at  the  points  (^1->2»n»Cj(_i)  • 


F'. ' /  =  R.H.S.  of  (3) F.  +  R.H.S.  of  (2,2)F 

I  ,K 


(3.24) 


Equation  ( 2 a3 ) F :  Since  there  are  now  solutions  available  at  two  previous 
C-positions,  a  3-point  derivative  in  c  is  used.  The  2-point  ^-derivative 
is  retained.  This  equation  is  used  at  the  points  (£-_i  .TifCj^)  and  has 
the  form 


F||k  -  R-H-S.  of  (2)F.  -  2(l-Ck)2(|F!ik-2F:>k.1+lF:>k.2)/ACk 


. 


, 


' 


. 
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25k<l-'k>KF;.k+1)(!  Fnk-2Fi,k-i+f  Fi,k-2> 


Fi;k  (|  Fi,k‘2Fi,k-l+  ?  Fi,k-2^/Ack  • 


(3.25) 


Equation  (3,3)F:  This  equation  is  used  over  the  interior  of  the  (£,?) 
plane.  Used  at  points  (£. >2 5  n  t  eniPl°ys  3-point  derivatives  in 
both  £  and  C- 


F!”  =  R.H.S.  of  (3) F.  +  R.H.S.  of  (2 ,3) F. 


(3.26) 


3.3c  The  Approximated  Energy  Equation 

The  energy  equation  also  requires  six  forms  for  use  at 
various  points  in  the  (£»?)  plane.  They  are 

Equation  ( 2) G :  For  points  (£._■]  ^C^q)  , 

G;;k  ■  -  (Fi.k^)6i.k-eitGi.ktFi.k-Fi-i.k)-tFi.k+i)(S.k-^-i.k)^i 

(3.27) 

Equation  (3)G:  For  points  (^-j >2 >n » ck<0^  * 


Gi;k  ■  -  (Fi.k+n)6,l.k-«l[Gi.k(|Fi,k-2Fl-l.k+?Fl-z.k) 


*Fi,k+1^2  Gi  ,k'2Gi-l  ,k+  2  Gi-2,k^/A^i 


(3.28) 


!. 
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Equation  (2,2)G:  For  points  (£.j_i  ,n»Ck=-| ) 


G:;k-R.H.S.  Of  (2)G.  =2(1-ck)2(G.jk-G.>k.iyACk 
"  2!k^‘?k^Fi  ,k+1^Gi  ,k'Gi  ,k-l^"Gi  ,k*Fi  ,k"Fi  ,k-l^/A?k 


(3.29) 


Equation  (3,2)G:  For  points  (5-j>2»rl»S|<=-| ) 


G!'k  =  R.H.S.  of  (3)G.  +  R.H.S.  of  (2,2)G 


(3.30) 


Equation  (2,3)G:  For  points  (£.=i 


G:;k  -  R.H.S.  of  (2)6.  =  2(l-Ck)2(|Gi>k-2Gi#k.1+^G.jk_2)/Ack 


-  2^k(! -Ck)  [(F^  #k+l )  (|  Gi,k"2Gi,k-l+  2  Gi  ,k-2^ 


Gi , k ^ 2  Fi,k”2Fi,k-l  +  2  Fi  ,k-2^/Ack  * 


(3.31) 


Equation  (3,3)G:  For  points  >2 »n »Ck>2) 


G!'k  =  R.H.S.  of  (3)G.  +  R.H.S.  of  (2,3)G. 


(3.32) 


3.4  Integration  in  the  and  Xr  Directions 

Computations  are  begun  by  finding  the  solution  to  the  initial 
steady-state  problem.  The  equations  to  be  solved  are  the  first  two 


■ 


‘ 
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in  section  3.3c.  Equation  (3.21)  is  used  for  i=l  and  equation  (3.22) 
for  i  ^2.  The  boundary  conditions  for  these  equations  are 


Fq  (Si  ,0)  =  0  , 
'i  'i 


F'  (5.  ,0)  =  jA-  Fi  ■(«,  ,0)-l  ;  5.  >  65.  , 

'i  'i  ^1.  'i  ‘i  'i  1 


F I  (^-1  »q-i  )  -  0  , 

'i  'i  ‘e 


F1  !^1 .  *nl  ^  0  * 

i  l  e 


(3.33) 


The  extra  condition  (on  F1')  is  required  because  the  value  of  ne  is 
a  function  of  £,  and  should  not  be  fixed  a  priori.  The  assumption  is 
also  made  that  free-stream  conditions  exist  at  the  leading  edge  of 
the  plate.  Thus 


F-,  (0,nq)  =  0  , 

i=0  1 

F-I  (O.n,)  =  0  ,  (3.34) 

i=0  1 

F-!'  (0,m)  =  0  . 

i=0  1 

As  a  consequence  it  is  not  necessary  to  compute  a  solution  at  £  =  0. 

The  equations  are  solved  by  marching  in  the  increasing  £- 
direction,  starting  with  £._i 


=0.1.  Six  steps  are  taken  with 


. 


8  II  I  gfi  I 
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A£..  =  -  5.. _i  =  0.1,  followed  by  seven  steps  with  .  =  0.2,  then 

six  steps  where  A£.j  =  0.4.  At  further  ^-positions  the  stepsize 
A£  =  0.8  is  used.  Marching  ceases  when  the  solution  obtained  is  within 
ten  percent  of  the  asymptotic  zero-slip  solution  obtained  by  Cohen  and 
Reshotko  [33].  This  point  is  reached  at  approximately  £  =  10. 

At  each  ^-position  the  unknown  boundary  condition  F!'(£.#0) 
must  be  found.  This  is  done  in  the  following  manner: 

(1)  A  guess  is  made  for  the  value  of  F!'(£..,0).  When  i=l,  the 
value  0.1  is  used.  For  i  >_  2,  the  value  found  at  the  (i-1) 
station  is  used  as  the  initial  guess. 

(2)  The  momentum  equation  is  integrated  over  the  q-variable 
until  a  value  of  n  is  reached  such  that  either  Fj(£..  ,n)  >_  0 
or  FI'(£. ,n)  <0-  If  integration  is  terminated  by  the  first 
condition,  the  guessed  value  of  F^'(C^.O)  was  too  large.  If 
termination  was  due  to  the  second  condition,  F^ ' ,0)  was 
guessed  too  small.  Based  on  this  information,  new  guesses 
are  made  at  intervals  of  AFV(^^,0)  =0.1. 

(3)  Once  the  true  value  of  Fl't^.  ,0)  has  been  bracketed  by  the 
above  method,  a  bisection  search  technique  is  used  to  obtain 
a  refined  estimate  of  Fl'f^.  ,0),  and  the  momentum  equation  is 
integrated  again  with  the  more  accurate  starting  value. 

(4)  The  search-integrate  procedure  is  repeated  until  the  conditions 
F! (^i  ,ne)  =  0  and  F! ' ,ne)  =  0  are  satisfied  to  four  places. 
This  requires  a  high  degree  of  accuracy  in  the  value  of 
F!'(S.,0)  due  to  error  propagation  in  the  integration  technique. 


-  n 

■ 

' 
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The  sensitivity  to  errors  is  inherent  in  the  Hartree-Womersly 
formulation. 

After  the  steady-state  momentum  equation  has  been  integrated 
out  to  £  =  10,  boundary  condition  (3.12)  is  used  to  obtain  the  velocity 
profiles  which  exist  immediately  after  acceleration  takes  place  (i.e. 
at  c  =  0).  Thus  we  calculate 


Uel  V2  U  ,  1/2  U  o  1/2 
F.  k_0U1,n.O)  =  (jjSi)  F.  ([rpl]  q.CT#]  n)  , 

l.KU  1  Ug3  I.  Ue3  1  Uel 

U  ,  U  ,  1/2  U  ,  1/2 

F!  k.0(5.-  .n,0)  =niiFl  £•  .  [rM]  n)  .  (3.35) 

1  ,K-U  1  Ug3  li  Ue3  1  U31 

U  i  3/2  U  I  1/2  U  3  1/2 
Fi:lk=o^i»n,°^  =  Fi ! ( ^i *  ^  * 


where  Ug3  is  arbitrarily  given  the  value  1.01  . 

We  now  have  the  solutions  which  are  assumed  to  exist  at 
C  =  0.  The  next  step  is  to  find  the  solutions  at  non-zero  values  of  c. 
The  method  used  is  similar  to  that  which  was  applied  for  non-zero  values 
of  £,  in  that  finite-difference  approximations  replace  the  ^-derivatives 
The  valid  equations  are  given  in  section  3.3c.  Equation  (3.23)  is 
used  for  the  case  when  (i=l,k=l),  and  equation  (3.24)  is  used  whenever 
we  have  ( i  >.  2 ,  k  =  1).  Expression  (3.25)  is  used  to  obtain  solutions 
at  the  points  (i=l,  k>2),  and  equation  (3.26)  represents  the  most 
general  case  where  (i  >_  2,  k  >.2).  Thus  equations  (3.23)  through  (3.26) 
are  used  to  march  slowly  in  the  increasing  c-direction.  The  stepsize 
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Ack  is  varied  as  follows:  A£k=1  =  0.05  and  Ack>2  =  0.1.  At  each  5- 
level  the  momentum  equation  is  integrated  over  0  <  £  £  10  using  the 
marching  technique  described  previously. 

The  distance  that  can  be  marched  in  the  timewise  direction 
is  limited  by  the  singularity  which  exists  at  t  =  0.5(i.e.  t=1).  Dis- 

w 

covered  by  Stewartson  [34]  and  later  confirmed  by  many  others  [35,36,1], 
the  singularity  causes  the  method  of  solution  to  become  unstable  in 
the  region  [18]  of 


c  =  - -  .  (3.36) 

2(F'+1 ) [(^y)  6+1]  +  (2-6) 

At  this  value  of  c  a  discontinuity  appears  [37]  in  the  f' distribution , 
and  the  present  method  becomes  incapable  of  proceeding  any  further  in 
the  c-di recti  on. 

Once  the  solution  to  the  momentum  equation  is  known  for  all 
0  <  £  <  10  and  0  <  c  <  L,  the  energy  equation  can  be  solved  with  the 
same  method. 

3.5  Integration  of  the  Ordinary  Differential  Equations 

Integration  of  the  O.D.E.'s  is  carried  out  by  means  of  a 
standard  integrating  subroutine  contained  in  IBM's  Scientific  Sub¬ 
routine  Package  [38].  This  subroutine,  called  DHPCG,  uses  Hamming's 
modified  predictor-corrector  method  with  double-precision  arithmetic 
to  integrate  a  system  of  first-order  O.D.E's  with  given  initial  values. 
It  is  therefore  necessary  to  replace  the  third-order  momentum  equation 
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with  three  coupled  first-order  equations,  and  the  second-order  energy 
equation  with  two  coupled  first-order  equations.  This  is  done  quite 
simply  by  defining  five  new  functions  such  that 


FI  =  F  , 


and 


F2 

F3 


dFl  _  c 
"  dn  '  F 

.  dF2  _  c 
dn 

G1  =  G  , 
dGl  _  r 

*r  G 


(3.37) 


(3.38) 


This  form  is  readily  accepted  by  subroutine  DHPCG.  Details  of  the 
procedure  can  be  found  in  Ref.  [38,39].  It  should  be  mentioned  here 
that  one  of  the  advantages  of  this  subroutine  is  that  it  will  maintain 
a  specified  accuracy  by  varying  the  stepsize  in  n- 

3.6  Errors  in  the  Smith-Cutter  Technique 

Since  the  ordinary  differential  equations  must  be  solved  by 
an  initial -value  procedure,  a  "shooting"  technique  is  used  in  order  to 
determine  the  initial  values  (unknown  wall  boundary  conditions)  F'*  and 
which  satisfy  the  known  free  stream  conditions.  The  correct  solution 
is  obtained  by  what  is  essentially  a  trial -and-error  method;  the  amount 
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of  error  introduced  is  dependent  on  the  degree  of  accuracy  with  which 
the  free-stream  conditions  are  satisfied.  In  the  present  case  an 
accuracy  of  four  places  was  specified.  Thus  the  maximum  error  intro¬ 
duced  by  the  shooting  technique  is  in  the  fifth  place  and  occurs  at  the 
edge  of  the  boundary  layer.  The  error  at  smaller  values  of  n  must  be 
less  since  about  ten-place  accuracy  in  the  wall  conditions  is  required 
to  obtain  four-place  accuracy  at  r\  . 

Errors  are  also  introduced  when  the  value  of  n  is  obtained. 

e 

As  the  edge  of  the  boundary  layer  is  approached,  subroutine  DHPCG  finds 
that  a  stepsize  of  An  =  0.1  is  small  enough  to  satisfy  the  error  re¬ 
quirements  imposed  on  the  integration  procedure.  Thus  the  value  of  ne 
is  accurate  only  to  one  decimal  place.  However,  the  results  show  that 
changing  the  value  of  ne  by  the  amount  An  affects  the  solutions  in  the 
fifth  place  only. 

Inaccuracies  due  to  round-off  errors  can  sometimes  be  significant. 
Double  precision  computations  using  the  very  fine  stepsize  An  =  0.01 
were  compared  with  solutions  obtained  using  the  stepsize  An  =  0.1;  agree¬ 
ment  existed  in  the  tenth  place. 

Smith  and  Clutter  found  that  the  boundary  layer  equations  as 
formulated  herein  would  become  extremely  sensitive  to  the  wall  conditions 
if  the  ratio  In/Aln  exceeded  25.  While  it  is  possible  to  control  this 
parameter  in  the  present  case,  there  is  an  additional  complication  which 
arises  as  a  result  of  the  unsteady  terms  8F'/3C  and  8G/9C-  These  terms 


have  the  form 
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whereas  all  other  finite-difference  terms  in  the  equations  can  be 
written  as 


Al 


n_  (  ) 

T  \  •  •  •  /  • 


n 


Even  if  is  kept  small  (unity,  for  example)  it  is  quite  possible 

for  1/ACj^  to  be  one  hundred  times  as  large.  The  result  is  that  the 
unsteady  term  dominates  the  differential  equation  in  such  a  way  that 
an  incorrect  value  of  F^'  or  produces  errors  which  increase  exponentially 
[27]  as  the  edge  of  the  boundary  layer  is  approached.  While  various 
numerical  methods  [40,41,42,43]  have  been  developed  to  eliminate  this 
problem,  a  less  exotic  procedure  known  as  the  Extended  Trajectory  of 
Integration  (ETI)  technique  is  used.  The  problem  is  not  removed,  it  is 
merely  circumvented  by  breaking  the  integration  into  a  number  of  shorter 
steps  over  which  the  errors  have  insufficient  time  to  get  out  of  hand. 

This  method  is  described  very  well  in  Ref.  [31], 

Smith  and  Clutter  [44]  use  a  surprisingly  simple  but  never¬ 
theless  excellent  method  to  investigate  the  question  of  instability  of 
solution  in  the  streamwise  direction.  They  do  this  by  solving  a 
problem  which  has  a  similar  solution.  A  similar  flow  is  one  of  the 
few  for  which  accurate  solutions  exist.  Because  the  flow  is  similar, 
the  solution  should  be  identical  at  all  stations  and  error  growth  is 
just  the  change  in  solutions. 
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In  the  present  case,  similarity  is  obtained  by  employing 
the  no-slip  condition  at  the  wall.  Then  the  steady-state  solution  be¬ 
comes  identical  with  that  obtained  by  Cohen  and  Reshotko  [33].  When 
both  types  of  solutions  are  compared,  agreement  to  four  places  is  ob¬ 
tained,  even  at  large  values  of  £.  Figure  2a  shows  errors  as  a 
function  of  £  at  the  n  =  0.1  level.  Figure  2b  shows  error  growth 
across  the  boundary  layer  at  the  station  £  =  4.4.  These  curves  are 
typical  and  show  that  the  errors  resulting  from  Smith  and  Clutter's 
technique  are  negligible  for  all  practical  purposes.  Since  no  similar 
solutions  exist  in  the  timewise  direction,  it  is  impossible  to  do  a 
similar  check  on  instability  along  the  £-axis.  We  can  only  hope  that 
the  accuracy  there  is  as  good. 

It  is  felt  that  the  solutions  presented  herein  are  accurate 
to  at  least  four  places.  This  upper  bound  on  the  accuracy  is  reached 
at  the  edge  of  the  boundary  layer;  accuracy  should  improve  at  lower 
values  of  n- 
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CHAPTER  IV 

APPLICATION  OF  THE  SOLUTIONS 

4.1  Shear  Stress  at  the  Wall 

The  shear  stress  at  the  wall  has  the  usual  form 

,  _  dUi 

^y=0 

Application  of  the  transformations  of  Chapter  II  yields 


(4.2) 


The  local  skin  friction  coefficient  may  be  defined  as 
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f 


(4.3) 


Combining  (4.2)  and  (4.3)  yields 


=  2 


(4.4) 


or 


(4.5) 
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We  now  introduce  a  Knudsen  number  based  on  fluid  properties 
evaluated  at  the  wall , 
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v  /u 
w'  w 


which  may  be  written  as 


(4.6) 


(4.7) 


Combining  this  relation  with  equation  (4.5),  the  result  is 


or 
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4.2  Heat  Transfer  at  the  Wall 

Aroesty  [13]  notes  that  the  energy  transferred  to  a  wall  from 
a  slightly  rarefied  gas  in  motion  is  given  by 
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where  q  ,  the  Fourier  heat  conduction  term,  has  the  usual  form 
w 


Q  =  k  Ml 
qw  Kw  9yly=0  • 


(4.11) 


and  u  $  ,  the  slip  work  term,  is  found  from 

w  w 


9u 


U  $  =  U  U  7 — 

w  w  Mw  w  9y'y=0 


(4.12) 


Since  Pr  =  CpP/k,  we  have 


C  y 

k  =  JOL 

Kw  Pr 


(4.13) 


Then 
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Pr  3y'y=0  w  w  y=0 


(4.14) 


This  equation  is  transformed  as  was  equation  (4.1).  The  result  is 


Q,.,  =  H 


'w 


P ,  a  (m+1)  V  1/2  C  e  . 
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(4.15) 


We  assume  that  Pr  =  1  and  note  that  if  the  stagnation  values  are  used 
as  the  reference  conditions. 
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Cp6ref  Cp0O  He 


(4.16) 


and 
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a  X  =  U  . 
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(4.17) 


Then  (4.15)  may  be  written  as 


Q„,  =  H  u 
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We  recall  that  equation  (2.45)  was 
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For  large  Mach  numbers  an  approximate  form  is 
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H 


=  2  . 


(4.19) 


Combining  (4.18)  and  (4.19)  we  have 


Pw  ae  (m+D  V  1/2 
Qw  =  He  U«  Pref  aref  [2vref  X  ]  Sw 


(4.20) 


Choosing  the  stagnation  values  for  reference,  and  using  expressions 
(A. 3),  (A. 5),  and  (A. 18),  equation  (4.20)  becomes 
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(4.21) 
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Defining  Nusselt  number  as 


Nu  = 


*S/0ref”0w^ 


we  obtain,  after  substituting  for  Qw  from  (4.21), 
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The  numerical  results  for  the  transient  contributions  to 
the  shear  stress  and  heat  transfer  at  the  wall  are  presented  in 
Appendix  D. 
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CHAPTER  V 

RESULTS  AND  CONCLUSIONS 
5.1  Discussion  of  the  Results 

The  transient  behavior  of  the  boundary  layer  on  a  flat  plate 
following  an  impulsive  acceleration  of  the  surrounding  fluid  has  been 
analysed.  Some  of  the  effects  of  low  density  have  been  included  by 
using  boundary  conditions  which  admit  a  slip  velocity  and  a  temperature 
jump  at  the  wall.  The  equations  governing  the  boundary  layer,  after 
being  rewritten  for  a  coordinate  system  suggested  by  the  boundary 
conditions,  have  been  solved  by  a  numerical  finite-difference  technique 
developed  by  Smith  and  Clutter  [27].  In  all  cases  the  unsteady  dis¬ 
turbances  were  initiated  by  an  impulsive  1%  increase  in  the  free  stream 
velocity. 

Several  special  cases  of  the  momentum  equation  (2.49)  have 
been  examined  by  other  workers.  Berard  [25]  studied  the  problem  of 
weak  interaction  (3=0)  steady-state  slip  flow.  Figures  3a  and  3b 
show  some  of  his  experimental  data  superimposed  on  theoretical  curves 
from  the  present  work.  The  best  agreement  was  obtained  when  curves 
corresponding  to  a  tangential  momentum  accommodation  coefficient 
a  =  0.6  were  compared  with  Berard 's  data  with  pressure  corrections. 

This  rather  low  value  for  a  corresponds  to  reflection  from  an  aluminum 
plate  covered  with  a  very  thin  film  of  oil.  The  oil  film,  which  was 
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observed  by  Berard,  is  more  "specular"  than  the  metal  surface  and 
hence  lowers  the  value  of  the  accommodation  coefficient. 

Aroesty  [13]  worked  on  the  problem  of  strong  interaction 
steady-state  slip  flow.  He  found  that  an  approximate  solution  for 
the  wall  shear  in  adiabatic  slip  flow  was 

V  =  f-'  -  eg  (5.i) 

W 

or,  in  terms  of  the  present  theory, 

f,','  =  0.765  -  |  .  (5.2) 

W  t, 

In  Figure  3c  this  solution  is  compared  with  that  from  the  present 
analysis.  Agreement  is  seen  to  be  good  except  in  the  region  of 
small  £,  where  Aroesty's  approximate  solution  predicts  a  separated 
fl  ow. 

It  is  seen  that  the  present  steady-state  theory  produces 
solutions  for  both  the  strong  and  weak  interaction  problems  which  com¬ 
pare  well  with  those  in  the  literature.  This  lends  support  to  the 
various  assumptions  made  during  the  course  of  the  analysis,  and  pro¬ 
vides  a  strong  base  from  which  the  assault  on  the  unsteady  problem 
can  be  launched. 

Examination  of  the  curves  in  Appendix  D  shows  that  the  un¬ 
steady  problem  has  been  solved  for  only  a  few  steps  in  either  the  £- 
or  £-  directions.  The  decision  to  limit  the  computed  solutions  was 
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based  on  purely  economic  reasons.  About  three  minutes  of  computer 
time  were  required  to  obtain  the  solution  at  each  point  in  the  £-£ 
plane.  The  original  intention  to  solve  at  25  £-steps  (to  £  =  10) 
and  10  c-steps  (to  £=0.5)  would  have  required  well  over  12  hours  of 
non-stop  computation.  It  is  felt  that,  since  the  method  of  solution 
has  been  developed,  the  solutions  at  larger  values  of  £  and  £  can  be 
obtained  with  relatively  little  difficulty. 

The  transient  contributions  to  the  velocity  and  shear 
functions  f'  and  f 1 1  are  plotted  in  Figures  4,  5a  and  5b.  The  tran¬ 
sients  are  observed  to  increase  with  £,  a  result  which  is  in  keeping 
with  the  fact  that  the  leading-edge  profile  does  not  change  with 
time.  Since  the  transients  are  referred  to  the  steady  slip  solu¬ 
tion,  the  ^-dependence  of  the  transients  is  not  related  to  the  £- 
dependence  of  the  slip  velocity,  but  is  instead  a  leading-edge 
effect  which  should  disappear  as  £  becomes  large.  Thus  there  will  be 
some  point  on  the  plate  beyond  which  the  transient  contribution  to 
the  velocity  and  shear  will  cease  to  be  a  function  of  £,  with  the 
result  that  the  transients  far  downstream  (£>10),  should  they  be  com¬ 
puted,  will  perhaps  approach  the  results  of  Rodkiewicz  and  Reshotko  [1]. 
This  trend  is  evidenced  by  the  increasing  curvature  of  the  slight 
hump  in  the  Af'  curves  at  low  values  of  n  as  £  is  increased. 

Figures  6a  and  6b  serve  to  illustrate  the  time-dependence 
of  the  drag  function.  The  drag  function  is  seen  to  decrease  mono- 
tonically  with  both  £  and  £,  varying  directly  with  the  square  of  the 
slip  velocity.  The  behavior  of  the  drag  near  the  leading  edge  is 
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open  to  question,  since  it  depends  so  much  on  the  assumed  shape  of 
the  leading-edge  profile.  However,  at  the  leading  edge  the  drag 
coefficient  defined  by  equation  (4.3)  can,  by  use  of  the  slip 
condition,  be  reduced  to 
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vr 


(5.3) 


which  is  the  value  corresponding  to  free  molecular  flow  [9].  Further¬ 
more,  examination  of  the  expression 


o  199  1/2 

cf  *  Cfw[sw +  }  >  (5-4) 

obtained  from  equations  (4.9)  and  (B.19),  shows  that  for  large  Mach 

numbers  the  factor  in  curly  brackets  is  greater  than  unity  over  at 

least  a  small  part  of  the  range  of  £.  For  large  £,  however,  the  low 

value  of  f'  forces  the  factor  to  become  less  than  one.  Qualitatively 
w 

this  behavior  agrees  precisely  with  that  described  by  Kogan  [8]  where i 
he  states  that  at  large  Mach  numbers,  as  the  Knudsen  number  decreases 
(increasing  £) ,  the  drag  and  the  heat  transfer  first  increase  above 
their  free-molecule  values,  reach  a  maximum,  and  then  decrease  to 
values  corresponding  to  a  continuous  medium.  A  detailed  study  of 
this  behavior  will  require  that  solutions  be  obtained  over  the  entire 
range  of  E,. 

The  transient  contributions  to  the  total  enthalpy  function 
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and  the  total  enthalpy  gradient  are  presented  in  Figures  7  and  8. 

The  transient  behavior  is  qualitatively  the  same  as  was  observed 
for  the  velocity  and  shear  functions,  with  the  addition  of  a  slight 
dependence  on  the  temperature  of  the  plate.  The  indications  are  that 
for  a  cold  plate  the  heat  transfer  in  the  interior  of  the  unsteady 
slipping  boundary  layer  is  less  sensitive  to  the  passage  of  time  than 
is  the  interior  region  of  the  boundary  layer  on  a  hot  plate.  Con¬ 
versely,  the  upper  region  of  the  cold  boundary  layer  seems  to  react 
more  quickly  to  the  passage  of  time  than  does  the  upper  region  of  a 
hot  boundary  layer.  This  behavior  contrasts  with  that  observed  by 
Gupta  and  Rodkiewicz  [18]  in  which  they  found  that  throughout  its 
thickness  a  cold  non-slipping  boundary  layer  reacted  more  slowly  than 
would  a  hot  zero-slip  boundary  layer.  They  attributed  this  effect  to 
the  fact  that,  for  the  hot  wall,  the  boundary-layer  density  is  less 
than  that  for  the  cold  wall,  rendering  the  hot-wall  boundary  layer  more 
susceptible  to  thermodynamic  changes  than  the  cold-wall  boundary  layer. 
Their  curves,  obtained  for  the  case  of  strong  interactions,  show  that 
the  sensitivity  difference  between  the  hot  and  the  cold  boundary  layer 
decreases  with  3-  A  linear  extrapolation  of  their  data  indicates  that 
for  the  case  of  weak  interactions  (3=0)  a  cold  boundary  layer  should 
be  more  sensitive  (have  larger  temporal  gradients)  than  a  hot  boundary 
layer.  This  agrees  with  the  present  findings,  at  least  in  the  outer 
regions  of  the  boundary  layer.  One  possible  explanation  for  the  in¬ 
sensitivity  at  the  interior  of  a  cold  slipping  boundary  layer  is  that 
the  gas  in  the  neighborhood  of  the  wall  (e.g.  0.1  <  n  <0.5)  is  com- 
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paratively  dense  so  that  its  temperature  is  insensitive  to  small 
fluctuations  in  the  heat  transfer  across  the  region.  The  absolute 
viscosity  and  the  heat  conduction  coefficient,  both  being  direct 
functions  of  temperature,  are  low  and  are  relatively  time  independent. 
Thus  the  heat  transfer  across  the  interior  region  remains  fairly 
constant,  and  most  of  the  transient  readjustments  take  place  in  the 
upper  levels  of  the  boundary  layer.  The  situation  immediately 
adjacent  to  the  wall  (n  <  0.1)  is  different.  There  the  gas  temperature 
is  low  enough  for  the  unsteady  variation  to  be  a  significant  fraction 
of  the  steady  state  value.  The  result  is  that  the  heat  transfer  at 
the  surface  of  a  cold  plate  is  more  sensitive  to  time  than  is  the  heat 
transfer  at  the  surface  of  a  hot  plate.  This  sensitivity  shows  up  as 
a  more  rapid  approach  to  the  steady-state  value,  as  illustrated  in 
Figures  9a  and  9b. 

The  steady-state  heat  transfer  distribution,  plotted  in 
Figure  9c,  is  seen  to  exhibit  a  monotonic  decrease  with  increasing  £, 
but  shows  no  sign  of  the  behavior  predicted  by  Kogan  [8].  This 
indicates  that  the  classical  slip  boundary  conditions  can  not  be 
relied  upon  to  provide  accurate  drag  and  heat  transfer  estimates  close 
to  the  leading  edge.  Therefore  the  drag  profile  given  in  Figure  6b 
should  be  used  with  caution,  as  the  matching  with  free-stream  conditions 
at  the  leading  edge  is  merely  a  coincidence,  and  is  not  a  result  of  the 
accuracy  of  the  boundary  conditions. 

Since  the  primary  purpose  of  the  present  work  was  to  explore 
the  transient  behavior  of  a  slipping  boundary  layer,  it  was  not  felt 
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necessary  to  include  plots  of  the  steady-state  velocity  and  temperature 
profiles.  Also,  solutions  to  the  steady  slip  problem  already  exist 
in  the  literature.  However,  for  reasons  of  completeness  this  and 
other  data  have  been  included  in  tabular  form,  and  may  be  found 
in  Appendix  E. 

Reproductions  of  the  computer  programs  used  to  solve  the 
energy  and  momentum  equations  are  presented  in  Appendix  F  with  the 
hope  that  they  may  be  of  some  assistance  to  the  reader. 

5.2  Concluding  Remarks 

The  present  results  may  be  used  to  obtain  the  history  of 
the  drag  and  heat  transfer  at  the  surface  of  a  flat  plate  in  slipping 
high-speed  flow.  Although  the  data  presented  is  quite  limited,  the 
computer  programs  are  such  that  very  few  changes  are  required  in  order 
to  produce  solutions  which  extend  almost  to  the  zero-slip  region. 

Once  the  solutions  for  the  entire  range  of  £  have  been 
completed  it  would  be  desirable  to  obtain  the  solution  to  cases  in  which 
there  are  strong  interactions.  This  would  of  course  require  that  the 
pressure  gradient  parameter  3  be  a  known  function  of  x  which  satisfies 
the  conditions  3 (small  £)  =  (y-l)/y  and  3(large  £)  =  0.  Similar  in¬ 
vestigations  for  other  airfoil  profiles  such  as  a  wedge,  circular 
cone  etc.  should  also  be  of  interest. 
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APPENDIX  A 

TRANSFORMATION  OF  THE  RAREFIED  GAS  BOUNDARY  CONDITIONS 
A.l  Velocity  at  the  Wall 

The  usual  boundary  condition  on  the  velocity  at  the  wall  of 
a  rarefied  gas  flowing  at  high  speed  over  a  surface  is  known  as  the 
"slip  velocity  boundary  condition".  Reference  [13]  shows  it  may  be 
written  as 


u(x,0,t)  =  AAW 


(A.l) 


where  A,  a  constant  for  a  given  gas-surface  combination,  is  computed 
from 


(A. 2) 


and  A  ,  a  function  of  x  and  t,  is  the  mean  free  path  of  a  gas  molecule 
adjacent  to  the  surface.  It  is  given  by 


(A. 3) 


If  the  transformations  of  equation  (2.33),  Chapter  II,  are 


applied  to  equation  (A.l),  then  the  transformed  boundary  condition  is 
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A2 


u(X,0,T)  =  AA 


w 


w  pref  aref 


3u(X,0,T) 

3Y 


(A. 4) 


When  we  define 


L  i  L(X.T) 


Vwae 


,i,refrjrefaref 


(A. 5) 


we  may  write 

u(X,0 ,T)  =  AArefL  3U(X’°’T>  .  (A. 6) 

At  first  glance  it  might  appear  that  L  depends  on  0W  and  thus 
creates  coupling  between  the  velocity  and  temperature  at  the  wall. 
Appendix  B  shows  why  this  need  not  be  so. 

We  introduce  the  incomplete  coordinate  transformations 


£  =  some  unknown  function  of  X  and  T  , 


n  =  [ 


(m+1 ) V  1/2 
-1 


2v  ,X 
ref 


(A. 7) 

(A. 8) 

(A.9) 


We  also  make  use  of  the  new  stream  function  f  defined  by 


2v  fX  V  1/2 

^  =  [  [ST]  1 


(A. 10) 
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Now 


u(X,Y,T) 


3^ 

3Y 


(A. 11) 


and 


If  =  Ve  f '(C.ri.T)  , 


(A. 12) 


or 


|f  =  aref  Me  f'U.n.t)  .  (A.13) 

Therefore 


u(X,Y,T)  =  aeMef'U,n,T)  .  (A. 14) 


Differentiating  with  respect  to  Y, 

3u  (m+1>  Ve 

W '  ae  Me  ^~2v  7"X^  f"U>n>T)  .  (A. 15) 

ref 


Substituting  (A. 14)  and  (A. 15)  into  (A. 6),  the  boundary  condition 


(m+1)  Vp  1/2 

aeMef'(5,0,T)  =  AArefLaeMe  [  ^ f"(C,0,x)  ,(A.16) 


becomes 


h 


A4 


which  reduces  to 


(m+l)  Vc  1/2 

f'(e,0,T)  =  AXrefL  [2—  ye]  f"U,0,x)  .  (A. 17) 


We  note  that  if  a  semi-similar  solution  was  desired,  wherein  f  =  f(n,i), 
the  above  slip  boundary  condition  would  be  functionally  inconsistent 
because  of  the  X-  and  T-dependence  on  the  right  hand  side. 

Since  the  independent  variable  £  is  still  unspecified, 
functional  consistency  will  be  obtained  if  we  let 


S  = 


VefL 


*w x 

L(m+1)  v; 


1/2 


(A. 18) 


Appendix  B  shows  that  this  definition  of  £  produces  an  independent 
variable  that  depends  only  on  X. 

From  (A. 18),  the  boundary  condition  becomes 


f'U.O.T)  =  |  f"(S,0,T)  . 


(A. 19) 


A. 2  Temperature  at  the  Wall 

Reference  [5]  gives  us  the  usual  form  for  the  boundary 
condition  on  the  gas  temperature  at  a  surface  immersed  in  high  speed 
slightly  rarefied  flow.  It  is 

e(x,0,t)  =  6p  +  , 


(A. 20) 


1  > 


A5 


where  0p  is  the  plate  temperature,  which  is  constant,  and  B  is  a  constant 
given  by 


b  -  (fir)  (p|r) 


(A. 21) 


When  equation  (A. 20)  is  transformed  as  was  (A.l),  the  result 


is 


6(X,0,T)  =  8  +  BA  ^e_  ieiMJl. 

P  w  Prefaref  3Y 


(A. 22) 


Introducing  L, 

9(X,0 ,T)  =  0p  +  BArefL  .  (A. 23) 

Dividing  through  by  0q  and  applying  (2.71),  we  obtain 

Q 

S(5,0,T)-f,2(£.0,T)  =  J-  +  BArefL  -^tS(C,0,T)-f,2(C.0,T)]  .  (A. 24) 

We  now  make  use  of  the  coordinate  transformations  found  in  section  A.l. 
Thus 


2  eD  (m+l )  V.  1/2 

S(5.0,t)  -  f^C.O.x)  =  +  BArefL  - jp]  [S'(C,0,t) 

(A. 25) 


-  2f '  (£  ,0  ,t)  f '  ' (£ ,0  ,x) ]  . 
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Substituting  (A. 18)  and  (A. 19)  into  the  above  equation,  then  rearranging, 
yields  the  boundary  condition  representing  a  temperature  jump  at  the 
wall , 


S(5,0,t)  =  +  |  S'(C,0,t)  +  iA-~gAP.L  f'  ,2(5,0,t) 

eo  ?  £ 

.  (A. 26) 

A. 3  Condition  for  an  Adiabatic  Wall 


From  Chapter  IV  we  know  that 

Q  =  H  p  a  /—  I  S'  U,0,x)  . 

Ww  eMw  w  ny  £  ' 

(A. 27) 

For  an  insulated  wall,  we  must  have  Qw  =  0,  and  therefore 


S'(5.0,T)ad  =  0  . 

(A. 28) 

(A. 28) 
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APPENDIX  B 

EVALUATION  OF  THE  QUANTITY  L, 

THE  PARAMETER  £,  AND  THE  KNUDSEN  NUMBER 

B.l  An  Examination  of  L 

We  wish  to  examine  the  quantity  L  more  closely.  Since 


(B.l) 


and 


'Vefprefaref 


(B.2) 


we  have 


or 


(B.3) 


Hayes  and  Probstein  [45]  note  that  at  high  temperatures 


B2 


1/2 


li  a  0 


(B.4) 


Thus 


1/2 


ew  ^ 


/ 


(B.5) 


or,  using  stagnation  conditions  as  the  reference, 


(B  .6) 


We  note  here  that,  although  (B.4)  is  inconsistent  with  equation  (2.20), 
its  use  is  very  convenient  in  that  it  renders  L  independent  of  the  gas 
temperature  at  the  wall.  Thus  L  no  longer  depends  on  time,  but  is  a 
function  of  X  only.  Since  both  (B.4)  and  (2.20)  merely  approximate 
the  true  dependence  of  viscosity  on  temperature,  the  inconsistency  is 
not  serious. 


If  the  Mach  number  is  large  an  approximate  form  for  (B.6) 


is 


L  = 


(B.  7) 
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B.2  An  Examination  of  £ 

We  wish  to  determine  how  £  depends  on  X.  The  definition  of 
equation  (2.34),  suggests  the  form 


K  = 


a  X1 


a 


(B  .8) 


Then 


a  = 


lref 


r2vref  X  "2  .-a 
L  L(m+1)  V  J  X 


(B  .9) 


But 


aref^e 


(B. 10) 


and  from  (B.7)  and  (B.10) 


Therefore 


(B.ll) 


[-2^4]1/2  x"“ 

Vef  aref  L  (m+1)  eXm 


a 


B4 
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1 

Vefaref 


1/2 

vref 


1/2  i  ■y'  m~7  (m'l)-“ 

e 


Since  we  know  that 


(B. 1 2) 


a  /  a(X) 


(B.13) 


we  must  have 


m 


-  2  (m_l )  -  a  =  0 


Using  the  relation 


(B.14) 


(B.15) 


we  have 


a  = 


1 

2-3  ‘ 


(B.16) 


Thus  £  has  the  form 


1 

2-3 

C  =  a  X  .  (B. 1 7) 


This  simple  expression  is  a  direct  result  of  assumption  (B.4),  which 
had  the  effect  of  removing  the  time-dependence  from  L,  and  thus  from 


\ 


Equation  (B.17)  can  be  differentiated  to  yield 


dl_  J_  £ 

dX  2-8  X  ' 


(B. 18) 


B.3  The  Knudsen  Number 

From  Equation  (4.6)  the  Knudsen  number  is  defined  as 


K  =  _ w_ 

n  Vuw  ' 


The  ratio  v  /u  can  be  thought  of  as  a  local  "characteristic  length" 
w  w 

for  slip  flow,  and  is  used  solely  for  convenience. 

The  grouping  A  /(v, Ju  )  has  a  rather  ubiquitous  nature.  As 
w  w  w 

written,  it  has  the  form  of  a  Knudsen  number,  the  significant  para¬ 
meter  for  slip  flow  [6].  However,  it  can  be  rearranged  to  u  A  ,/y  . 

w  w  w 

a  form  of  Reynolds  number,  or,  using  equation  (B.l),  to  /  tty/ 2'  M  , 

w 

a  direct  function  of  the  local  Mach  number  at  the  wall. 

An  expression  for  the  Knudsen  number  can  be  obtained  in  terms 
of  the  numerical  solution  as  follows: 
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Kn  = 


u. 


/Tjx1  W 

2  aw 
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or 
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f  ' 
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CO  + 


III 
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2)S  -  111 

i  pw  2 


M  2f'2] 
e  w  J 


-1/2 


(B. 19) 


At  the  leading  edge  this  becomes 


e  * 


(B.20) 


whereas  further  downstream  an  approximate  form  is 


w 


(S  -f'2) 
x  w  w  ' 


1/2 


(B.21 ) 


which  is  independent  of  the  free-stream  conditions  so  long  as  the  Mach 
number  is  large.  Essentially  equation  (B.19)  describes  how  the  slip 
effects  decay  as  the  no-slip  region  is  approached. 
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APPENDIX  C 

CORRELATION  WITH  AROESTY'S  WORK 


We  want  to  learn  the  relationship  between  Aroesty's  perturba¬ 
tion  parameter  e  and  the  modified  x  coordinate  We  have  [13] 


£ 


\.PWU« 


/IT 


(c.l) 


and 


€  = 


J _ 

VefL 


rfVref  X 


1/2 


(C.2) 


We  note  that 


where 


n^(n  as  used  by  Aroesty)  = 


P  dy  , 


(C.3) 


C  = 


x 
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pw  ^  U6  dx 


(C.4) 


and 


(C.5) 


C2 


In  this  analysis 


n  =  [ 


(m+1)  Vn  1/2  a, 
-] 


2vref  X 


ref 


ref 


dy  , 


(C.6) 


where 


f  a  P0 

x  =  cmm  (p-5-) dx 

J  aref  ref 


(C.7) 


Thus 


JL-  r (m+1)  V/2  ae  1  1  ^ 

nA  l2vrefX  aref  pref  Ue 


(C.8) 


Using  (C.l)  in  (C.8),  we  obtain 
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Employing  equation  (C.2),  the  result  is 


_n_  _  _]_ 
nA 


Returning  attention  to  (C.8),  we  see  that 


(m+1)  V  1/2  a 

_n  -  r  e1  e  _J _ L 
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where  (C.4)  and  (C.5)  have  been  used  with  the  relations 
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and  Reference  [46]  gives  us 
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■ 


Therefore 
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(C.20) 


But  from  (C.16)  we  have 
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Equating  the  right-hand  side  of  (C.20)  with  the  right-hand  side  of 
(C.21),  and  substituting  the  result  into  equation  (C.18),  we  obtain 
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If  C  *  constant  (an  approximation  made  by  Aroesty)  then  we 


can  write 
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Evaluating  the  integral  in  the  numerator  produces 


JL 

nA 


3y_1  m(l-  ^-  +  1 

(m+l)(l-  jfcl)  x  4y  4y 


1/2 


3v,  {m+1 )  ( 1  -  ^4~) 

(nt+1 )  (1-  x _ 

_[  3y  1 

(m+1)  (1-  ^jj~)  x  4y 


1/2 


] 


=  1  . 

We  therefore  have  the  simple  relations 

n  =  nA  * 


(C.24) 


(C. 25) 


and 
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APPENDIX  D 

CURVES  FOR  CHAPTERS  III  AND  IV 

-containing  a  graphical  presentation  of  error  growth,  compari¬ 
son  with  the  work  of  others,  and  the  numerical  solutions. 
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FIGURE  2b.  Error  Growth 


D4 


o 


a> 

CD 


FIGURE  3a.  Gas  Velocity  at  the  Surface 

a  Semi-Infinite  Flat  F 
in  Steady  Slip  Flow 
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FIGURE  3b.  Velocity  Profile  on  a  Se  m  i  -  I  n  f  i  n  i  t  e 

Flat  Plate  in  Steady  Slip  Flow 
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FIGURE  3c.  Wall  Shear  Along  a  Semi  -  Infinite  Flat  Plate 

in  Steady  Adiabatic  Slip  Flow 
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FIGURE  4.  Transient  Contribution  to  the  Velocity  Functi 


D8 


o 


D  9a 


A(£fK 


D  9b 


FIGURE  6a.  Transient  Contribution 
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FIGURE  6b.  Steady  -  State 

the  Drag 
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Enthalpy  Function  at  £  =  0.3  with 
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FIGURE  9b.  Transient  Contribution  to  the  Heat  Transfer  Function  with 
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APPENDIX  E 

TABLES  OF  STEADY  STATE  SOLUTIONS 
-the  steady-state  profiles  of  f',  f ' 1 ,  S,  and  S' 


form. 


in  tabular 


f=  0.1  ,  A*?  =  o.l 


STEADY  STATE  (WITH  SLIP)  VELOCITY  DISTRIBUTION: 


0.92364 

0.93253 

0.94071 

0.94820 

0.95501 

0.96117 

O.96669 

0.97161 

0.97596 

0.97979 

0.98312 

O.986OO 

0.98847 

0.99058 

0.99235 

0.99384 

0.99508 

0.99610 

0.99693 

0.99760 

0.99814 

0.99857 

0.99891 

0.99918 

0.99938 

0.99954 

0.99966 

0.99975 

0.99982 

0.99987 

0.99991 

0.99993 

0.99995 

0.99997 

0.99998 

0.99998 

0.99999 

0.99999 

0.99999 

0.99999 

0.99999 

1.00000 

STEADY 

STATE  (WITH 

SLIP)  SHEAR 

DISTRIBUTION 

* 

0.09236 

0.08532 

0.07834 

0.07148 

0.06479 

0.05834 

0.05216 

0.04631 

0.04082 

0.03571 

0.03100 

0.02670 

0.02282 

0.01935 

0.01627 

0.01357 

0.01122 

0.00920 

0.00748 

0.00603 

0.00482 

0.00382 

0.00300 

0.00233 

0.00180 

0.00137 

0.00104 

0.00078 

0.00058 

0.00043 

0.00031 

0.00023 

0.00016 

0.00011 

0.00008 

0.00006 

0.00004 

0.00003 

0.00002 

0.00001 

0.00001 

0.00000 

2  ,  A 17  =  c 

i.l 

STEADY 

STATE  (WITH 

SLIP)  VELOCITY  DISTRIBUTION: 

0.85260 

0.86907 

0.88435 

0.89845 

0.91136 

0.92310 

0.93372 

0.94324 

0.95171 

0.95920 

0.96577 

0.97148 

0.97642 

0.98064 

0.98422 

0.98724 

0.98976 

0.99185 

0.99356 

0.99495 

0.99608 

0.99697 

0.99768 

0.99824 

0.99868 

0.99901 

0.99927 

0.99946 

0.99961 

0.99972 

0.99980 

0.99986 

0.99990 

0.99993 

0.99995 

0.99996 

0.99997 

0.99998 

0.99999 

0.99999 

0.99999 

0.99999 

0.99999 

1.00000 

STEADY 

STATE  (WITH 

SLIP)  SHEAR 

DISTRIBUTION 

• 

• 

0.17052 

O.I5876 

0.14688 

0.13500 

0.12325 

0.11174 

0.10059 

0.08989 

0.07973 

0.07018 

0.06129 

0.053H 

0.04565 

0.03891 

0.03290 

0.02758 

0.02293 

O.OI889 

0.01544 

0.01250 

0.01004 

0.00799 

0.00630 

0.00492 

O.OO381 

0.00293 

0.00223 

0.00168 

0.00125 

0.00093 

0.00068 

0.00049 

0.00035 

0.00025 

0.00018 

0.00012 

0.00009 

0.00006 

0.00004 

0.00003 

0.00002 

0.00001 

0.00001 

0.00000 

• 

• 

- 

. 

• 

• 

, 

f  0.5  ,  AU  =  0. 1 


E3 


STEADY 

STATE  (WITH 

SLIP)  VELOCITY  DISTRIBUTION: 

0.78719 

0.81006 

0.83141 

0.85122 

0.86949 

0.88622 

0.90142 

0.91515 

0.92746 

0.93840 

0.94805 

0.95651 

O.96385 

0.97017 

0.97557 

0.98015 

0.98399 

0.98719 

0.98983 

0.99199 

0.99375 

0.99515 

0.99628 

0.99716 

0.99786 

0.99839 

0.99881 

0.99912 

0.99936 

0.99953 

0.99966 

0.99976 

0.99983 

0.99988 

0.99992 

0.9999^ 

0.99996 

0.99997 

0.99998 

0.99998 

0.99999 

0.99999 

0.99999 

0.99999 

0.99999 

1 . 00000 

STEADY 

STATE  (WITH 

SLIP)  SHEAR 

DISTRIBUTION 

• 

0.23616 

0.22117 

0.20586 

0.19040 

0.17493 

0.19561 

0.14461 

0.13006 

0.11610 

0.10285 

0.09040 

0.07882 

0.06817 

0.05847 

0.04973 

0.04194 

0.03507 

0.02907 

0.02388 

0.01945 

0.01570 

0.01256 

0.00996 

O.OO783 

0.00609 

0.00470 

0.00359 

0.00272 

0.00204 

0.00152 

0.00112 

0.00082 

0.00059 

0.00042 

0.00030 

0.00021 

0.00015 

0.00010 

0.00007 

0.00005 

0.00003 

0.00002 

0.00001 

0.00001 

0.00001 

0.00000 

1  ,  At?  =  0 

'•i  .  ep/ea 

=  1 

STEADY 

STATE  (WITH 

SLIP)  ENTHALPY  DISTRIBUTION: 

0.92473 

0.93349 

0.94156 

0.94894 

0.95565 

0.96172 

0.96716 

0.97201 

0.97631 

0.98007 

0.98336 

0.98620 

0.98864 

0.98071 

0.99246 

0.99393 

0.99515 

0.99616 

0.99698 

0.99764 

0.99817 

O.9986O 

0.99893 

0.99919 

0.99940 

0.99955 

0.99967 

0.99976 

0.99983 

0.99988 

0.99991 

0.99994 

0.99996 

0.99997 

0.99998 

0.99999 

0.99999 

0.99999 

1.00000 

1.00000 

1.00000 

1.00000 

STEADY  STATE  (WITH  SLIP)  S-PRIME  DISTRIBUTION: 


0.09105 

0.08411 

0.07722 

0.07046 

O.O6387 

0.05751 

0.05142 

0.04565 

0.04024 

0.03520 

0.03056 

0.02632 

0.02250 

0.01907 

0.01604 

0.01338 

0.01106 

0.00907 

0.00738 

0.00595 

0.00475 

0.00376 

0.00295 

0.00230 

0.00177 

O.OOI36 

0.00103 

0.00077 

0.00057 

0.00042 

0.00031 

0.00022 

0.00016 

0.00011 

0.00008 

0.00006 

0.00004 

0.00003 

0.00002 

0.00001 

0.00001 

0.00001 

. 

■ 

• 

• 

f  =  0.2  ,  A»7  =  0.1  , 


@p/@o  ~  1 


E4 


S 


s 


STEADY  STATE  (WITH  SLIP)  ENTHALPY  D I STR I  BUT  ION : 


0.86861 

0.88257 

O.89568 

0.90790 

0.91921 

0.92959 

0.93904 

0.94759 

0.95525 

0.96206 

0.96806 

0.97331 

0.97787 

0.98178 

0.98512 

0.98794 

0.99030 

0.99226 

0.99388 

0.99519 

0.99626 

0.997H 

0.99779 

0.99832 

0.99874 

0.99906 

0.99930 

0.99949 

0.99963 

0.99973 

0.99981 

0.99986 

0.99991 

0.99993 

0.99996 

O.QQ997 

0.99998 

0.99999 

0.99999 

0.99999 

1 . 00000 

1.00000 

1.00000 

1.00000 

STEADY 

STATE  (WITH 

SLIP)  S-PRIME  DISTRIBUTION: 

0.14363 

0.13541 

0.12671 

0.11767 

0.10844 

0.09916 

0.08996 

0.08097 

0.07228 

0.06400 

0.05619 

0.04893 

0.04224 

0.03616 

0.03068 

0.02581 

0.02153 

0.01779 

0.01458 

O.OII83 

0.00952 

0.00759 

0.00600 

0.00470 

0.00364 

c. 00280 

0.00213 

0.00161 

0.00120 

O.OOO89 

c. 00065 

0.00048 

0.00034 

0.00024 

0.00017 

0.00012 

0.00008 

0.00006 

0.00004 

0.00003 

0.00002 

0.00001 

0.00001 

0.00001 

3  ,  a  17=  0.1  ,  ep/e, 

0=  1 

STEADY 

STATE  (WITH 

SLIP)  ENTHALPY  DISTRIBUTION: 

0.83054 

0.84707 

O.86283 

0.87773 

0.89172 

0.90475 

0.91677 

0.92779 

0.93779 

0.94679 

0.95483 

0.96194 

0.96818 

0.97360 

0.97827 

0.98225 

O.98562 

0.98845 

0.99079 

0.99272 

0.99430 

0.99557 

0.99658 

0.99739 

0.99802 

0.99852 

0.99889 

0.99918 

0.99940 

0.99957 

0.99969 

0.99978 

0.99984 

0.99989 

0.99993 

0.99995 

0.99997 

0.99998 

0.99999 

0.99999 

0.99999 

1.00000 

1.00000 

1.00000 

1.00000 

STEADY 

STATE  (WITH 

;  slip)  s- 

-PRIME  DISTRIBUTION: 

0.16888 

0.16157 

0.15341 

0.14455 

0.13513 

0.12529 

0.11522 

0.10507 

0.09500 

0.08516 

0.07567 

0 . 06664 

0.05818 

0.05033 

0.04315 

0.03666 

0.03086 

0.02574 

0.02128 

0.01742 

0.01413 

0.01136 

0.00905 

0.00714 

0.00558 

0.00432 

0.00331 

0.00252 

0.00189 

0.00141 

0.00104 

0.00076 

0.00055 

0.00040 

0.00028 

0.00020 

0.00014 

0.00010 

0.00007 

0.00004 

0.00003 

0.00002 

0.00001 

0.00001 

0.00001 

. 


* 

f  =  o.i  ,  A/?  =  o.i  , 


0p/0o  =  0.5 


E5 


STEADY  STATE  (WITH  SLIP)  ENTHALPY  DISTRIBUTION: 


0.89165 

0.90426 

0.91587 

0.92649 

0.93616 

0.94489 

0.95273 

0.95971 

O.96589 

0.97132 

0.97604 

0.98013 

0.98364 

0.98663 

0.98915 

0.99127 

0.99302 

0.99447 

0.99565 

0.99660 

0.99757 

0.99798 

0.99846 

0.99884 

0.99913 

0.99935 

0.99952 

0.99965 

0.99975 

0.99982 

0.99987 

0.99991 

0.99994 

0.99996 

0.99997 

0.99998 

0.99999 

0.99999 

1.00000 

1.00000 

1.00000 

1 . 00000 

STEADY 

STATE  (WITH 

SLIP)  S-PRIME  DISTRIBUTION: 

O.I3IO7 

0.12108 

0.11117 

0.10143 

0.09194 

0.08278 

0.07402 

0.06572 

0.05792 

0.05067 

0.04399 

0.03789 

0.03239 

0.02746 

0.02309 

0.01926 

0.01 592 

0.01306 

0.01062 

0.00856 

0.00684 

0.00542 

0.00425 

0.00331 

0.00255 

0.00195 

0.00148 

0.00111 

O.OOO83 

0.00061 

0.00044 

0.00032 

0.00023 

0.00016 

0.00012 

0.00008 

0.00006 

0.00004 

0.00003 

0.00002 

0.00001 

0.00001 

f=o.2  ,  aj?=  o.i  ,  0P/eo  =  °* 5 


STEADY  STATE  (WITH  SLIP)  ENTHALPY  DISTRIBUTION: 


0.80429 

0.82546 

0.84525 

0.86364 

0.88058 

O.89609 

0.91018 

0.92288 

0.93423 

0.94431 

0.95317 

0.96091 

0.96761 

0.97337 

0.97826 

0.98240 

O.98585 

0.98872 

0.99108 

0.99300 

0.99455 

0.99580 

0.99678 

0.99756 

0.99816 

0.99863 

0.99899 

0.99926 

0.99946 

0.99961 

0.99972 

O.9998O 

0.99986 

0.99991 

0.99994 

0.99996 

0.99997 

0.99998 

0.99999 

0.99999 

1.00000 

1 . 00000 

1.00000 

1.00000 

STEADY 

STATE  (WITH 

SLIP)  S-PRIME  DISTRIBUTION: 

0.21832 

0.20489 

0.19094 

O.17667 

0.16228 

0.14794 

0.13385 

0.12017 

0.10705 

0.09457 

0.08289 

0.07205 

0.06211 

0.05309 

0.04500 

0.03781 

0.03150 

0.02601 

0.02129 

0.01727 

O.OI388 

0.01106 

0.00874 

0.00684 

0.00530 

0.00407 

c. 00310 

c . 00234 

0.00175 

0.00129 

0.00095 

0.00069 

0.00050 

0.00035 

0.00025 

0.00018 

0.00012 

0.00008 

0.00006 

0.00004 

0.00003 

0.00002 

0.00001 

0.00001 

- 

• 

\ 

f  =  o.j  ,  Aj?=o.i, 


E6 


0P/0o  =  °-5 


STEADY  STATE  (WITH  SLIP)  ENTHALPY  DISTRIBUTION: 


0.73697 

O.76361 

0.78880 

0.81245 

0.83449 

0.85488 

0.87359 

0.89064 

0.90605 

0.91985 

C. 93212 

0.94294 

0.95239 

0.96058 

O.9676I 

0.97359 

0.97864 

0.98286 

0.98636 

0.98924 

0.99157 

0.993^6 

0.99496 

0.99616 

0.99709 

0.99782 

0.99838 

0.99880 

0.99913 

0.99937 

0.99955 

0.99968 

0.99977 

0.99984 

0.99989 

0.99993 

0.99995 

0.99997 

0.99998 

0.99999 

0.99999 

1.00000 

1.00000 

1.00000 

1.00000 

STEADY  STATE  (WITH  SLIP)  S-PRIME  DISTRIBUTION: 


0.27339 

0.25932 

0.24430 

0.22852 

0.21219 

0.19554 

0.17881 

0.16221 

0.14597 

0.13027 

0.11529 

0.10117 

0.08802 

0.07591 

0.06490 

0.05499 

S  0.02098 

0.04618 

0.03844 

0.03170 

0.02591 

0.01684 

0.01339 

O.OIO55 

O.OO823 

0.00637 

0.00488 

0.00370 

0.00278 

0.00207 

0.00153 

0.00112 

0.00081 

O.OOO58 

0.00041 

0.00029 

0.00020 

0.00014 

0.00010 

0.00006 

0.00004 

0.00003 

0.00002 

0.00001 

0.00001 

f=  0.1  ,  A  7?  =  0 

.1  , 

=  0 

0 

STEADY 

STATE  (WITH 

SLIP)  ENTHALPY  DISTRIBUTION: 

O.85856 

0.87502 

0.89018 

0.90405 

O.91667 

0.92807 

0 . 9  385c 

0.94741 

0.95548 

O.96256 

0.96875 

0.97407 

O.97865 

0.98255 

0.98584 

o.9p86o 

O.99089 

0.99278 

0.99432 

O.99556 

S  0.99657 

0.99736 

0.99799 

0.99848 

O.99885 

0.99916 

0.99938 

0.99955 

0.99967 

0.99977 

0.99983 

O.99988 

O.99992 

O.99994 

0.99996 

0.99998 

0.99998 

O.99999 

0.99999 

1 . ooooc 

1 . 00000 

1 . boooo 

STEADY 

STATE 

(WITH 

SLIP)  S-PRIME  DISTRIBUTION: 

0.17109 

0.15805 

0.14511 

0.15240 

0.12002 

0.10806 

0.09662 

0.08579 

0.07561 

0.06614 

0.05742 

0.04946 

0.04227 

0.03584 

0.03014 

0.02513 

0.02079 

0.01705 

0.01386 

0.01117 

r .00893 

0. 00707 

0.00555 

0.00432 

0.00333 

0.00255 

0.00195 

0.00145 

0.00108 

0.00079 

0.00058 

0. 00042 

0.00030 

0.00021 

0.00015 

0.00011 

0.00007 

0.00005 

0.00003 

0.00002 

0.00C02 

0.00001 

-- 

£=0.2,  A  17  =  0. 1  , 


E  7 


S 


S 


/ 


£  =  o. 


S 


S 


/ 


STEADY  STATE  (WITH  SLIP)  ENTHALPY  DISTRIBUTION: 


0.73997 

O.76835 

0.79483 

0.81937 

0.84195 

0.86260 

0.88132 

O.89817 

0.91322 

O.92656 

0.93828 

0.94851 

0.95736 

0.96495 

0.97141 

0.97685 

0.98141 

O.98518 

O.98828 

0.99081 

0.99285 

0.99449 

0.99578 

0.99680 

0.99759 

0.99820 

0.99867 

0.99903 

0.99929 

0.99949 

0.9996^ 

0.99974 

0.99982 

0.99988 

0.99992 

0.99994 

0.99996 

0.99998 

0.99998 

0.99999 

0.99999 

1.00000 

1 . 00000 

1.00000 

STEADY 

STATE  (WITH 

SLIP)  S-PRIME  DISTRIBUTION: 

0.29301 

0.27437 

0.25518 

0.23567 

0.21611 

0.19672 

0.17774 

0.15937 

0.14179 

0.12515 

0.10959 

0.09518 

0.08198 

0.07003 

0.05931 

0.04981 

0.04146 

0.03422 

0.02800 

0.02270 

0.01825 

0.01453 

0.01147 

O.OO898 

O.OO696 

0.00534 

0.00407 

0.00307 

0.00229 

0.00170 

0.00124 

0.00090 

O.OOO65 

0.00046 

0.00033 

0.00023 

0.00016 

0.00011 

0.00007 

0.00005 

0.00003 

0.00002 

0.00001 

0.00001 

3  1  A»7  =  0 

.1  ,  eje 

0=  0 

STEADY 

STATE  (WITH 

SLIP)  ENTHALPY  DISTRIBUTION: 

0.64339 

0.68015 

0.71477 

0.74716 

0.77725 

0.80501 

0.83041 

0.85349 

0.87430 

O.8929I 

0.90941 

0.92393 

O.9366O 

0.94755 

O.95695 

0.96495 

O.97166 

O.97728 

0.98193 

0.98575 

O.98885 

0.99135 

0.99335 

0.99492 

0.99615 

0.99712 

O.99786 

0.9984 2 

O.99885 

0.99917 

0.99940 

0.99957 

0.99970 

O.99979 

O.99986 

O.99990 

0.99993 

O.99996 

O.99997 

O.99998 

0.99999 

O.99999 

1 . 00000 

1.00000 

1.00000 

STEADY 

STATE  (WITH 

SLIP)  S-PRIME  DISTRIBUTION: 

0.37790 

O.35708 

0.33518 

0.31248 

0.28925 

0.26579 

0.24239 

0.21935 

0.19693 

0.17538 

0.15491 

O.13569 

O.II786 

0.10149 

0.08664 

0.07332 

0.06150 

0.05113 

0.04213 

0.03440 

0.02783 

0.02232 

0.01773 

0.01396 

0.01089 

0.00841 

0.00644 

0.00489 

O.OO367 

0.00273 

0.00202 

0.00147 

0.00107 

O.OOO76 

0.00054 

0.00038 

0.00027 

0.00018 

0.00013 

0.00008 

0.00006 

0.00004 

0.00002 

0.00002 

0.00001 

. 

. 

. 

\ 

FI 


APPENDIX  F 

THE  COMPUTER  PROGRAMS 

-  written  in  FORTRAN  IV  for  an  IBM  360/67  computer  using  the 
MTS  system.  Data  generated  by  the  programs  is  stored  on  magnetic 
disks  using  the  MTS  FILE  system. 


c 

c_ 


MAIN  PROGRAM 


C 


c 

c 


c 


c 

c. 


c 


c 


c 


c 

c 

c 


c 


IMMICIT  RL.ALA8  (A-H.O-Z) 

PI  ML  NS  I  UN .  F(  3.  /I  ,3)  ,  FN  (  3, 71 . 3  )  .TNN(  3,  7  L,  3L,_E  NT1N  (  7 1.) . 

Di MFNSI UN  PHMT ( j ) , FUNC ( 3 ) « DERF 11), AUX ( 16,3) 

DIMENSION  STfJRb<  3,71  ,2,3)  ,  STORE  A  (  26,71  ,3  )  ,  STURfcB  (3,71  ,  3  ) 
DIME  NS  ION  FNNS( J),  I  SOL ( 2 ) 

COMMON  F ,FN,FNN,FNNN,  I  SOL , NE /E CT D A 1  /LQN,X.DX, T , D  T , G 
T1CMMUN  /INTQAT/  IX  ,MPS _ 

EXTERNAL  FCT,OUTP 


I  0-0 

10=  l 


F.1  The  Momentum  Equation 
Program 

NS=2  _ _ . _ __________ _ 

NS=3 

UE3= 1,01  00 


KT  -  2 
K  T=  3 
KT- 1 
KT- O 


I  LX-b 
I  EX  =  3 


1BSR=0 
NO I M=-3 

PRMT(  2)  =700 

PRMT  C  3  )  .=  (>,  U1Q _ 

PRMT (4) =50-2 

1 TERS=30 
I  TER 5=4  0 
ITERS =30 

IF  (  KT,  GT  ,0  ) _ IT  E  H_S  ~3  0 

DO  6  L-  1 , 3 

DO  o  J-  1  ,  7\ 

DO  6  1=1, 26 _ _ 

STORE  A (  I  , J,L)=ODO 
IF  (I  ,1)1*3)  GO  TO  S 
ST URL (  I  , J ,  1  , L  J-ODG 
STORE (  I , J ,2 , L ) =0  DO 
STOREb! I , J , L ) =ODU _ 


S  CONTINUE 


' 


■ 


' 

■ 


<•>  CONTINUE 

C: 

DO  10  J-l ,71 

c  _ 

F  NNN (  J  )  ~0D 0 
C 

DU  10  1-1,3 

C 

DO  10  K  =  1  *  3 
C 

FNN (  I »  J  ,  K ) =0  00 
FN(  I  »  J  ,  K ) — ODO 
F (  I  ,  J  ,K)  =  QDG 
C. 

10  CONTINUE 
C 


IF  (KT.NF.1 )  GO  TO  12 
C 


READ  (1)  ST OR FA 


f" 

1  5 

CONTI NUE 

DU  11  1—1,3 

DO  11  L-  1  »  3 

DO  11  J=1 , 71 

c 

STORE (  I  •  J  *  1  ,L 1-5 TORE AC  1 , J , L ) 

c 

1  1 

CCNTI NUF 

IF  (NS.EQ.l)  GO  TO  22 

c 

12 

CONTINUE 

c 

c 

IF  (KT  *LT ,2 , UR , KT,GT,3 )  GO  TO  14 

READ  (1)  STOREA 

IF  (KT.E0.2)  READ  (2)  STORE 8 

IF  (KT.EQ.3)  READ  (3)  STOREB 

c 

DO  13  1-1,3 

DO  13  L  = 1  «  3 

DO  1J  J=l,71 

c 

f" 

STORE!  I  ,  J,  2, L) -STOREA!  I  , J,L) 
STORE!  I  , J  ,1 »  L ) - ST  ORE 0 (  I  , J,L) 

1  3 

CONTINUE 

c 

14 

CUNT  I NUE 

c 

DO  20  1-1,3 

c 

FNNS! I ) =0D0 

c 

20 

CONTINUE 

C 


F3 


5 


. 


c 


c_ 

c. 

c 


G-bOO/ IDO 
G=1 .  400 
TMACCO-O .600 
_ TMACCU-  1QQ_ 

DX-0.1D0 
DN=0. IDO 
DT-0.0 1 2500 
DT=0 .0500 


FNN W=0  .  IDO 
DFNNW  =  0.  IDO 


£ 


C 


C 


22 


C 


C. 

C 


c 

c 


c 


PI  —  3.  141  59  2  6b  35  8  ->  7  9  JD  0 
DLX=DX 

OTEL-Or 
X  — ODO 
T-ODO 

Cl-  ( 2  DO  —  T MACCO ) /TMACCO 
NE=71 


Df=8D0*DT  EE 

IF  (KT.GE.12)  DT  —  0.01 DO 
IF  (KT.LE.9)  DT=4D0*DTE£ 
IF  (KT.LE.5)  DT=2DO*OTEE 
IF  (KT.LE.2)  DT=DT  LE¬ 


IF  (KT.GE.12)  T-G.bDO+DFLOAT( KT- I  1  ) *DT 
IF  (KT.LE.li)  T -0 • 3D0  +  DFLG AT ( KT— 9 ) *DT 
IF  (KT.LE.9)  T  =  0.  1DO  +  DFLOAT ( KT— 5 )  *DT 
IF  (KT.LL.5)  T=0. 1 DO+DFLOAT (KT— 2 ) #DT 

IF  (KT.LE.2)  T-DFLQAT(KT) *DT 


CONTINUE 
DO  180  IX- 1,1  LX 

PRMT ( 1 ) =000 

IF  (IX.EQ.26)  PRMT ( 4 ) =5D— 5 

IF  ( KT.EQ.O . AND.NS.NE. 1 )  GO  TO  26 

DO  23  K  =  2  *  3 

KKK=K— 1 
L=K 

IF  (NS.EQ.l)  L=  1 

IF  ( NS.EQ.l )  KKK=1 

DO  23  J=l .71 

F ( l • J «L)=STUPE( IX , J.KKK, 1 ) 

FN(  1  . J.L)=STORE(  IX.  J.KKK.2) _ 

FNN (  1 ,J.L)=STOPL<  IX,  J.KKK. 3) 


23  CONTINUE 


. 


* 


F5 


c 

26 

£ 

C 


C 


£ 

C 


c 


c 


c 

'30 

c 

'35 

C 


C 


C 


C 

37 

C 


C 


C 


IF  (NS.EQ.t)  GO  TO  98 
CONT I NUE 


IF  (IX.EQ.l)  GO  TO  35 


I  J=1 


IF  (  IX .EQ.2 .OR. I X.E0.7. OR .  IX .EQ. 1 4 .UR.  IX  .EQ. 20 >  IJ  =  2 


DO  30  I  1=  I  J  »  2 

1=4-1  I 
I JI =1-1 


DO  30  J  =  1  .71 

F(I *J*1) =F (IJI.J.l) 

FN ( I , J  *  1  ) -FN (  I J I . J ,  1  ) 

FNN{  I  *  J,  1  )=FNN(  I JI* J,  1  > 

CONTINUE 

CONT I NUE 

MDPTCR=0 

LF=0 

_ 

DX=8D0*DEX 

IF  { I X*LE .19)  DX=4D0*DEX 
IF  (IX.LE.13)  DX=2D0*DEX 
IF  (IX.LE.fe)  DX=DEX 
X=X+DX 

E QN=2D0 

IF  (  IX.GT. I .AND. KT .EQ.O)  EQN= 3D0 

IF  ( IX.EQ.l .AND. KT.EQ. 1 .OR. KT.EQ. 12)  EQN=22D0 

IF  < IX. GT. 1 . AND .KT .EU. 1 . 0R.K1  . EQ. 12 )  EGN  =  32D0 

IF  ( IX.EQ.l .AND. KT.GT. 1 )  EQN=23D0 _ 

IF  ( IX.GT . 1 . AND.KT.GT . 1 )  EQN=33D0 

CONTINUE 


I QSR=0 
I  SOLI  1 )=Q 


MPS=  1 

IF  (MDPTCR.EQ.l )  MPS=11 

IF  (KT. EQ.O. AND. IX.GT. 2)  FNNW=2D0 *F NN ( 2 .  1  *  1  ) -F NN ( 3 ,  I  *  1 ) 
IF  (KT.EQ.l)  FNNW=Q .99D0*FNN( 1  *  1*2) 

IF  (KT.GT.l)  FNNW=2D0*FNN< 1* 1*2)— FNNI1, 1,3) 

IF  (MPS.EQ.ll)  FNNW=FNN<  1  *  1  1,  1) 


C 


■ 


. 

' 


c 

c 

40 

c 


c 


c 

c 

c 

_ 

c 


c 


F6 

IF  (KT.EQ.O.AND.  IX.GT  .2 )  DFNNW=D AE3 S ( FNN ( 2 ,  I ,  1 ) -F NN ( 3 * l  • 1  )  ) / 1 00 0 
IF  (KT.EQ.l)  DFNNW=0  #0 1  DO  *FNN(  1*1*2) 

IF  (KT.GT.l)  DFNNW=DAB S(FNN ( 1 , 1 , 2) -FNN( 1 , t .3 ) ) /I ODO 
■IF  (  MPS.,  EQ«  11  )  DFNNW—  O  .O  I  QOtFNNM 

IF  (MPS.tQ.ll )  FNNA=FNNW 

CONTI NUt 

IF  (MDPTCR.EQ  .0  )  LF=LFtl 

IF  (MDPTCR. EQ . 1 )  LF1~LF1+1 

I  SOL ( 2 ) =0 

DERF ( 1 )=0. IDO 
DERF ( 2)=0.3D0 
PERF(3)=Q.6D0 

IF  (MDPTCR.EQ. 1 )  GO  TO  63 

FUNC(  3  ) -FNNW/ 0.900 

CONTINUE 

FUNC( 3)=0.9D0*FUNC(3 > 

FUNC( 2 )=C 1 *FUNC<  3 ) /X-1D0 
IF  (FUNC( 2) . GT .000 )  GO  TO  60 
FUNC ( 1 ) — ODO 


GO  TO  66 
C 

63  CONTINUE 
C 

PRMT( 1 ) — l DO 
FUNC ( 3 ) =FNNW 

FUNC( 2)=FN( 1*11*1) 

FUNC ( I )  =  F(  1,11,1) 

C 

66  CONTINUE 
C 

FNN  S  (  2  )=FUNC(  3)  

C 

CALL  DHPCG ( PRMT  *  FUNC  *  DERF  »NDIM,  I HLF  »  F  C  T  *  OUTP  *  AUX ) 

C 

IF  (IHLF.GE.il)  WRITE  (6,70)  I HLF 
70  FORMAT  (*  • /T 1 0 ,  *  I HLF= • , 1 2/ ) 

IF  (IHLF.GE.il)  STOP 

C 

ERR0R=5D-S 
E RROR=  50— 6 

IF  (IX. GT  •  24  )  E«ROR=5L»-4 

C 

_ I  F  (  DAO S ( FUNC( 2 )  )  .LT. ERROR.  AND.  DABS ( FUNC( 3 )  )  . L T . ERROR)  GO  TO  90 

IF  (LF1.EQ.50)  GO  TO  90 

IF  ( LF. EQ . 1  • OR. (LF 1 .EQ. 1  , AND. MPS.EQ . I  1  )  >  I  SOL (  l ) = I  SOL ( 2 ) 

IF  ( (LF.NE  .  I  .OR • (LF1 . NE . 1  . AND • MPS. EQ .  1 1 )  )  .AND. (  ISOL(  1  )  .NE. 


■ 


F7 

1 ISOLC2) •OR.I0SR.GT.O) )  GO  TO  80 
C 


FNNS(  I  SOL ( 2 )  )  =FNNS(2) 

_ IF _ l  I  SOI  (  ?  >  .  EQ  .  1  )  FNNW=FIMN«I  +  DFNN  W 

IF  ( I  SOL ( 2)  . EQ.3 )  FNNW=FN NW— OF N N W 

C 

GO  TO  40 

C 

c 

80 

CONTINUE 

I BSR= I 

FNNS (  I  SOL (2)  )  =FNNS ( 2 ) 

FNNW=(FNNS( 1 ) +FNNS ( 3) ) /2DG 

IF  (LF+LF1 .£Q . ITERS )  MDPTCR'-l 

IF  ( LF+LF 1 .EQ. I TERS)  GO  TO  37 

c 

GO  TO  40 

c 

90 

CONTINUE 

c 

IF  (NE.EQ.71>  GO  TO  98 

c 

JS=NE+ 1 

DO  9  5  J-JSF71 

c 

DEL=1D0— DFLGAT (7 1  - J ) /DFLOAT ( 71-NE) 

F  (  1  «  J  •  1  )=OEL*F(  i  •  J  •  1  1 

FN (  1  ,  J,  1 )  —  DEL  *FN (  1 , J,  1  ) 

FNN(  1  ,  J,  1 ) =DEL*FNN<  1  .  J  .  1  ) 

FNNN ( J ) =DEL*FNNN  (  J  ) 

c 

95 

CONTI NUE 

c 

96 

CONTI NUE 

c 

IF  (MPS.EQ.I)  GO  TO  98 

c 

DELFNN=(FNN(  1.11*1  )-FNNA)/2DQ 

c 

DO  97  J  =  2  * 2 0 

c 

ETA=DFLOAT( J-l )/ 1 ODO 

IF  (J.LE.ll)  FNN{  1 , J, 1  ) =FNN (  1  , J .  1 ) F LI A*DELFNN 

IF  (J.GT.ll)  FNN{ 1. J, 1)-FNN( 1* J. 1 ) - ( 2O0-LT A ) 4DELFNN 

c 

97 

CONTINUE 

c 

90 

CONTI NUE 

IF  (NS.EQ.2)  GO  TO  165 

c 

DO  100  J  =  .1  .  NE 

C 


C 


FNd.Jtl  )  =f"N  l  l.J,  1  )  +  lD0 


. 


. 


F8 

100  CONTINUE 

c 

FSE  TA=DN*DFLOAT  (  NE-  1  ) 

C 

IF  (IX.EQ.l)  WRITE  (6.110) 

110  FORMAT  ( ' 1 • ) 

C 

T  S—  T 

IF  ( KT. EQ. 0 . AND. NS .EQ.  3 )  T-10D0 

_ WRITE  (6.1201  X.FSETA.T _ _ _ _ _ 

T=TS 

120  FORMAT  (  * 6 '  .  T  1  0  .  '  X~*  ,F 4. 1  , T20,  •  ET A( EPEE  STREAM  )-•  »F3. 1  • 

1  T45,  *  T  =  *  , F 6 . 4 ) 

C 

WRI TE  (5.130) 

--  1  30 _ FORMAT  (  * _ '  //  /T  10.*  UN  S  T  E  AD  Y  STATE  (WITH  SLIP)  VELOCITY  PISTRIBU 1 Li 

1  «n  :  »  / ) 

C 

WRITE  (6,140)  (FN(  1,J»1).J=1. NE ) 

140  FORMAT  (•  •.T2.SF23.5) 

C 

WRITE  (6.  ISO)  

150  FORMAT  (•  ' ///T 1 0 .' UNSTEADY  STATE  (WITH  SLIP)  SHEAR  DISTRIBUTION:1 

1  /) 


C 


c 


156 


C 


c 

c 

c 

160 

1  65 
C 

c 


c 

c 

c 


WRITE  (6,140)  ( FNN( 1 , J, 1 ) , J= 1 ,NE ) 

LFS  =  LF+I F  1 


WRITE  (6,155)  LF.LFl.LFS 

FORMAT  (»  '/////T10,  'CONVERGED  IN  *  ,  12 , *  +  *  • I 2  •  *  =  *  ,  I  2,  *  ITERATIONS, 

1  •  ) 

WRITE  (6,110)  

DO  160  J=1,NE 
F  N (  1  , J ,  1 ) -F N (  1 , J ,  1  ) - 1  DO 

CONTI  NUE 

CONTINUE 

IF  (KT.LU.O)  GO  TO  171 


JI  =  2 

IF  (KT.EQ.Q.OR.KT.EQ. 2 . OR . K T . EQ . 5 .UH.KT.EQ.9 .QR.KT  ,EQ. 1 1  )  J I _ 

DO  170  K-l.JI 
DO  170  1 , 71 

STORE (  I  X  .  J , K ,  1 ) -E  (  1  , J ,K  ) _ 

STORE (  IX.J.K,2)“FN(  l.J.K) 

STORE (  IX, J.K,3)=FNN(  1  , J.K) 


C 


■■■ 


. 


1  70 


CONTINUE 


C 

1  71 

C 

C 

C. 


C 

1  72 

C 

1  76 
C 

1  80 
C 


c 

c 

c 

1  90 

c 

c 

c 


c 

2  00 
c 


c 

c 

210 

c 


c 

P.PO 


CONTINUE 


IF  (KT.GT.O)  GO  TO  176 

DO  172  J=l,71 

STORE A (  I  X , J,  1  )  =F (  1  *  J,  1  ) 
STORE  A (  IX, J,2)=FN(  1 * J,  1  ) 
STORE A (  IXiJ,3)  =  FNN(  1  *  J  *  1  ) 

CONT I NUE 

CONTI NUE 


CONTINUE 

IF  (IO.EQ.O)  GO  TO  220 

IF  (KT.NE.O)  GO  TO  200 

IF  (NS.EQ.l)  GO  TO  220 

IF  ( NS  «EQ.2)  GO  TO  190 

WRITE  (5)  STORFA 

GO  TO  220 

CONTI NUE 


CALL  STATE2IST  OR  t A ) 
WR ITF  (  1  )  STORE A 


NS-  1 

GO  TO  1 5 

CONTINUE 

DO  210  1—1*3 

DO.  21  Q.  L  =  1  t  3 

DO  2  10  J  =  1  *  7 l 

STOR  EB  {  I  ,  J  ,  L  )  =S  T  ORE  (  1  *  J  *  1  *  L) 
CONTINUE 


IF  (KT.EQ.l)  WRITE  (2)  STOREB 
IF  (KT.EQ.2)  WRITE  (3)  STORES 
IF  (KT.EQ.3)  WRITE  (A)  STORES 

CONTI NUE 


C 


STOP 

END 


•,  ■ 


c 


c. 


c 


c 


c 


c. 


c 

c 


c 

c 


c 


c 


c 

c 


c 


c 


SUBROUTINE  FCT(ETA,FUNC,OERF) 

I  MPL  I  C  I  T  RE  AL  (A-H.O-Z) 

-PI  MENS  IUN  F  (  3  .  ^ NMU  .7  1,  3l»f  NNM(  71  j 

DIMENSION  FT( 6 1  *  2)  *  ARG( 20) • VAL( 40  > 

DIMENSION  FUNC  (  3  )  •  DERF  (  3  ) 

COMMON  F  ,  FN , FNN, FNNN 

COMMON  /FCTD AT/tQN  »X  »  OX  *  T  *  DT  *  G 

COMMON  / INTDAT/I X.MPS 

ERR-5D-5 

ERR-5D-6 

IF  (IX.GT.24)  ERR=SD-5 

N1-MPSF60 

N2  =  20 
N3  =  6  1 
EPS=EHR 


BGN-ODO 

-LF (MB3 « Efl«  I  1  )  BGN=1DJ> 

DERF ( I ) =FUNC( 2 ) 

D£RF(2)  =  FUNC<  3) 


I NTERP=  1 

DO  3  J=MPS »  7 1 

N-J 

ET  AC  HK  =  0 • 1 D  0  +DFLOA  T  (  N-i  } 
PET  A- FT  A— FT  ACHK 


IF  (DABS(DETA) .LE.SD-10)  INTERP-0 
IF  (INTERP.EQ.O)  E  I  A—  ET  ACHK. 

IF  (INTERP.EQ.O)  GO  TO  6 
IF  ( DET  A.LE.ODO)  GO  TO  o 


3  CONTINUE 

b  CONTINUE 

F  Ml  X—  F  (  2  .  N.  1  ) 
FNM1X=FN( 2.N. 1 ) 

FM1T=F(  t  *  N • 2 ) 
FNM1 T=FN< 1 .N .2 ) 

FM2X=F ( 3 • N  «  1  ) 
FNM2X=F N( 3,N,  1  ) 


C 


F M2T=F< 1 ,N.3) 
FNM2T-FN ( 1 ,N ,3 ) 


c 


FI  1 


IF  ( INTERP.EQ.O)  GO  TO  50 
C 

DO  IQ  J-MPS.N1 

c 

J J- J-MPS+  i 
FT(JJ,1)=F(2,J,1) 

FT ( JJ ,2)=FN( 2, J,  1  ) 

c 

10  CONTINUE 

C 

CALL  CAT  5E ( ET  A , EGN,  0 . I  DO . FT • N3, Z , ARG ♦ VAL,N2) 

C 

CALL  DAHI ( ET A, AR G, V AL • P.N2.EPS, I LR ) 

C 


FMIX^P 

C 


20 


C 

C 


c 

30 

C 

c 


c 


c 


4  0 


c 


IF  (IER.GT.l)  WRITE  (6,20)  ILR 

FORMAT  {•  ,/T10f  '  IER=» , I  1  *  *  FOR  FM1X.V) 

IF  (IER.EQ.l)  EPS— 1 Q DO* EPS 
IF  (  IER.EQ.l)  GO  TO  10 
JiPS=ERR 

IF  (  IER.GT.0 )  STOP 

DO  30  J  -MPS  *  N  1 

J J=J-MPS+ 1 

FTC JJ ,  1  )=FN(2«J, 1  ) 

FTC  J J ,2) =FNN( 2  «  J  «  1  ) 

CONT I NUE 

CALL  DATSECETA ,BGN,  0.  1D0,F  T,N3,  2,  ARG*VAL»N2  ) 


CALL  DAHI (ETA, ARG, VAL ,P,N2,EPS, I ER J 


FNM1 X=P 


IF  (IER.GT.l)  WRITE  (6,40)  1 ER 

FORMAT  («  */T 10  ,  «  IER= •  , I  I ,  »  FOR  FNM1X .  «✓ ) 

IF  (IER.EQ.l)  EP  S—l ODO*£PS 
IF  (IER.EQ.l)  GO  TO  30 
EPS* ERR 

IF  (  IER  ,GT  .0  )  STOP 


50  CONTINUE 

C 

IF  (EQN.NE.200)  GO  TO  60 
C 

DERF  (  3)  — —  (  FUNG  (  1  )  +ET  A  )  *FiJNC(3  )  -X*  (  (FUNC(  1  )  -FM1X  )*FUNC  (3  ) 
1  — ( FUNC ( 2  ) — FNM IX ) *(FUNC< 2  )  +  lDO ) ) /DA 

C 


c 


RETURN 


60  CONTINUE 


B 


* 


FI  2 
c 

IF  ( INTERP.EU.0.CR.FCN.EQ.3DU)  GU  TO  110 
C 

DO  70  J=MPS.N1 

C 

J J= J-MPS+ 1 

FT<  JJ , 1 ) =F<  1  ,  J  ,  2  ) 

FT(JJ,2)=FN(  1  ,  J  *  2  ) 

C 

 7.Q CON  TIN UE 

C 

CALL  OAT  SE ( ET  A  , BGN  *  0 • 1  DO , FT , N 3 , 2 • ARG , VAL , N2 ) 

C 

CALL  D AH  I  (ETA  ,  ARG, VAL • P , N2 * EPS.  ILR) 

C 


c 


30 


c 

c 


c 

90 

c 

c 


FM1T-P 


IF  (IER.GT.l)  WRITE  (6,80)  IER 

FORMAT  (*  « /T10* •  IER=»  • I  1  •  •  FOR  EMIT.*/) 

IF  (IER.EQ.l)  EPS=10D0*EPS 
IF  (IER.EQ.l)  GO  TO  70 
EP.S.= ERH 

IF  (IER.GT.O)  STOP 

OO  90  J=MPS,Nl 

JJ-J-MPS+ 1 
-FI.ljJLJ-a.1  J^NtUJig) 

FT (  J  J  *  2  ) =FNN(  1  ,  J,2) 

CONTI NUE 

CALL  DAT  SE  (  ET  A  *  BGN  *  0  .  1  DO  *  F  T  *  N3 , 2.  ,  AR  G  ,  V  AL  •  N2  ) 


C 

c 


12)32 


c 

1  1  0 

c 

c 


CALL  C  AH I  (ETA,  ARG* VAL, P.N2.EPS.  IER) 
FNM1T=P 


IF  (IER.GT.l)  WRITE  (6,100)  IER 

_ f  Q.EM  AJ— .  1 .  * _ *  /T  10,*  IER  ~  *  ,  I  1  _ EAR  _FNML1,»ZI  _ _ _ 

IF  (IER.EQ.l)  EPS=10D0*£PS 
IF  (IER.EQ.l)  GO  TO  90 
EPS=ERR 

I F  (  IER  .GT . 0 )  STOP 

CONTI  NUF _ _ _ 

IF  (EQN.NE.22C0 )  GO  TO  120 

IF  (MPS.EQ.ll  .AND.N.LU.  1  1)  DERF  (  3)  ~FNN  N(  1  1  ) 

OERF (3 ) =- ( FUNC (  1  )  +ET  A  )  *F UNC ( 3 ) -X* (  (FUNC (  1)-FM1X) *FUNC (  3 ) 

1  —  ( F UNC ( 2  )  -F NM1X)* (FUNC ( 2  )  +•  1 00 ) ) / 1 > X _ 

2+200*  (  1  DO-T  )  **2*(  FUNC  (  2  ) — FNM1  T  )  /DT-2U0*  T ♦  (  1D0-T  ) 

3* (  (FUNC (2  )-FNMl T ) * (FUNC (2 )  +  1D0 )  —  ( FUNC (  1  ) -FM 1 T ) *FUNC ( 3 )  )/0T 


C 


RETURN 


FI  3 


C 

120 

JL 

c. 

c 


c 

1  30 

c 

c 


CONTINUE 

IF  ( INTERP.FQ.O )  GO  TO  170 

UC1  130  J=  MRS  *  N  1 

J J- J-MPS+ 1 
FTIJJ.1  )-Ff3.  .1.1  ) 

FT(JJ,2)=FN( 3, J, 1 ) 

CONTINUE 

CALL  DATSE(ETA ,BGN, 0.  1 D 0 , F T , N3 , 2 ♦ AR G, V AL , N2  ) 


C 

C 


140 


C 


c 


c 

150 

c 

c 

c 


CALL  D AH I (ETA.ARG, VAL , P , N2 , EPS , I ER ) 

FM2X=P 

IF  (IER.GT.l)  WRITE  (6,140)  I ER 

FORMAT  C  *  /  T  1 0  , 1  IE  R= *,11,*  FOR  FM2X  .  « /  ) 

IF  (IER.EQ.l)  EPS=10DQ*EPS 
IF  (IER.EQ.l)  GO  TO  130 
EPS=ERR 

IF  (IER.GT.O)  STOP 

DO  ISO  J=MPS.N1 

JJ= J-MPS+1 
F  T ( J  J  •  1  )=FN(  3, J,  1 ) 

FT ( JJ ,2>=FNN( 3, J, 1 ) 

CONTINUE 

CALL  DATSEtET A , 6GN, 0 • 1 D 0 , F T , N3 , 2 , AR G * V AL , N2 > 
CALL  DAHI (ETA, Aft G , V AL , P , N2 , EPS , I ER ) 

FNM2X=P 


C 


1  60 


C 

1  70 

C 

C 


IF  (IER.GT.l)  WRITE  (6,160)  I ER 

FORMAT  (•  • / T  1 0 ,  •  I ER= **11**  FOR  FNM2X •  * / ) 

IF  (IER.EQ.l)  EPS=10D0#EPS 

IF  (IER.EQ.l)  GO  TO  150 

EPS=ERR 

I F  (  IER.GT.O )  STOP 
CONTI NUE 

IF  (EQN.NE.300)  GO  TO  180 


DERF  (  3  )  (  FU  NC  (  1  )  FETA  )  *P  UNC  (  3)  -X*  (  (  1  .  SDO*FUNC(  I  )  -2D0*FM1  X  +  0 .500 

1  *FM2X ) 4FUNC ( 3  )- (  1  .5D0*FUNC(2 ) -2D 0 *FNM 1 X F 0 . 5D0 
2*FNM2X)<‘(FUNC(2)  +  100)  )/DX 


' 


, 


. 


* 


. 


c 


RETURN 


c 

_ LiiO _ CUNT  I  NUF  _ _ _ _ _ _ _ _ 

c 

If  ( EON. NE. 32  CO )  GO  TO  190 

c 

DERF  (  3  )=-  (FONC  (  1  >  +ET A )  *FUNC(  3  )-X*(  (  1  .500  *FUNC  (  1  ) -2D0*F  M  1  X  +  G  •  5D0 
1 *FM2X) *FUNC ( 3 ) - ( 1 •SD0*FUNC<2 ) -2D0*FNM 1 X+0.5D0 

_ 2_*-irNM2X)*<FUNC<2)-flL)0)  )/px _ 

3  +  200  *  (  1 DO-T ) 4*2* ( FUNC ( 2 ) — FNM 1 T ) /DT-2D0*T*(  1D0-T ) 

4* ( (FUNC( 2 ) -FNM 1 T )  *  (FUNC (2 )+l DO ) - (FUNC ( 1 J-FMl T ) *FU NC ( 3 ) )/DT 
C 

RETURN 

C 

190  CONTINUE 

C 


IF  ( I NTERP.  EQ • 0 )  GO  TO  240 
C 

DO  200  J=MPS,N1 
C 

J  J- J -MP3 + 1 

F  T(  J  J  ,  1  )  =F  (  1  ,  J  ,3  ) 

FT< J J  ,2  )=FN<  1  , J.3) 

C 

200  CONTINUE 
C 

_ CALL  D  AT  SE  (  F.T  A  ,  EGN  ,  O.  1  DO,  FT  ,  N3,  2  ,ARG,V  At  ,N  21 _ 

C 

CALL  DAH I  (ETA*  ARG *  V  AL  *  P  »  N2  » EPS*  I ER ) 

C 


C 


FM2T=P 


210 


C 


c 

22  0 
C 

C 

C 

c 


IF  (IER.GT.1)  WRITE  (6,210)  ICR 

FORMAT  (•  VT10,  *  IER-'  .  I  1  ,  •  FOR  FM2T  *  *  /  ) 

IF  (IER.EQ.l)  EPS=10D0*EPS 

IF  (IER.EQ.l)  GO  TO  200 

EPS=ERR 

IF  (IER.GT.O)  STOP 


DO  220  J-MPS.N1 

J J- J-MPS+ 1 

FT( JJ, 1 )=FN( 1 , J, 3) 

FT( J J. 2) =FNN( 1 , J, 3> 

CONTINUE 

CALL  CATSE( ET A, BGN, 0 . IDO, FT  ,N3.2  * ARG*  VAL  ,N2) 
CALL  DAH  I  ( ETA ,ARG, VAL ,  P, N2.EPS, I ER ) 

FNM2T  =P 


IF  (IER.GT.1)  WRITE  (6,230)  IER 


' 


. 


230  FORMAT  (•  • /T 10 • • 1EH“ •  • I  1 •  •  FOR  FNM2T •  •  /  ) 

IF  <  IER  .  EQ  .  1  )  EP5=10D0*EPS 
IF  (IER.EQ.l)  GO  TO  220 
EPS=ERR  

I F  (  IER  .GT  .0  )  STOP 


240  CONTINUE 
C 

IF  (fcQN.NE.23C0)  GCJ  TO  250 

jC  

DERF ( 3 )=- (FUNC(  1  )  +ET  A  )  *FUNC( 3 )-X*(  <FUNC(  1  )-FM IX) *F UNC ( 3) 

1  — ( F  UNC ( 2 ) -F  N  M 1 X ) * (FUNC ( 2 )  +  1  DO ) ) /DX 

2  +2 DO  * ( 1  DO  — T )**£*<  1 . 5D 0*FUNC ( 2 ) -2D0 *FNM 1 T+ 0 . 5D0*F NM2 T ) /DT 
3-2D0*T*< 1D0-T )*( ( 1 . 5DO*FUN C ( 2 ) —200 *FNM 1 T + 0 • 50 0 *F NM2 T ) 

4 * (FUNC (  2) + 1  DO ) -(  1 .5D0*FUNC (1  ) -2 DO *FM l T  +  0 .5D0*FM2T ) *FUNC< 3 )  ) /DT 
C 


RETURN 

C 

250  CONTINUE 

C 

DERF ( 3 ) =— ( FUNC ( 1 ) +ETA ) *FUNC ( 3 ) -X* ( ( 1. SDO*FUNC< 1 )-2D0*FMl X+0.5D0 
l  *FM2X)»FUNC(3)-(  I  .500  *FUNC  (  2  )-2D0*FNM  1  X+0.  SDQ 

2*FNM2X) * ( FUNC ( 2 ) + 1 DO ) )/UX 

3  +2 DO  4(1 DO-T ) **2*(  1 .5D0*FUNC< 2 )-2D0*FNM 1 T  +  0 . 5D0*FNM2T )/DT 
4-2D0*T* ( 1 DO-T ) * (  <  1 .SDC*FUNC<  2  >-2D0 *FNM I T  +  0 .5D0*FNM 2T ) 

5  4  (  F  UN  C  (  2  )  FIDO  )-(  1 .5D0*FUNC(I  ) -2D0*FM  IT-FO  .  5DQ*FM2T  }  *FUNC  <3  )  ) /OT 

c 

RETURN 

END 


SLBROUT I NE  GU TP ( ETA ,FUNC • DERF , I HLF,  ND I M *  PRMT  ) 
IMPLICIT  REAL*8  (A-H, O-Z) 


— nu  ii  ui  v-ii* — i — i  j  t  i  i  tJi  u  m  jj.iitj;  ir^  i  a  »  i.  i  a  ^  \  z  i.  j _ 

DIMENSION  PRMT ( 5 ) ,FUNC(3) ,OERF( 3) , IS0L(2 ) 

c 

COMMON  F ,FN , F  NN, FNNN* I  SOL , N 

COMMON  /INTDAT/ IX.MPS 

c 

NX=  3 

NX  =  4 

NX=  1 

NX=  2 

c 

INTERP=1 

c 

DO  3  J=1 ,71 

c 

n 

z 

c 

ET  ACHK=0 •  1 DO *D FLO AT ( N— 1  ) 

DET  A=ET A— ET  ACHK 

c 

IF  (OABS(DETA) .LE.5D-1  0)  INTERP=0 

IF  ( INTERP.EQ .0 )  ET  A  =  ET  ACHK 

IF  ( INTERP.EQ. 0 >  GO  TO  6 

IF  ( DETA  .LE .0D0)  GO  TO  6 

c 

3 

CONTINUE 

C 

6 

CONTINUE 

c 

IF  (  IX.GE.NX)  IG-l 

10=1 

10=0 

c 

IF  (  I NTERP.NE  .0  )  RETURN 

c 

c 

F( 1 ,N, 1 ) =FUNC ( 1 ) 

FN(  1 ,N, 1 ) =FUNC( 2  ) 

FNN( 1 , N, 1 )=FUNC( 3 ) 

FNNN( N) =DERF( 3) 

c 

IF  ( FUNC(  2 ) .GT.-5D-6 )  I  SOL  ( 2)  =3 

IF  (FUNC(3) ,LT,SD— 6)  ISOL(2)=l  _  _  . 

C 

c 


1  0 


1 


20 


IF  l I  SOLI  2)  .NE.O)  PRMT (5) -I  DO 
IF  (  ICUEQ.O  )  GO  TO  30 

WRITE  (6*10)  ETA,N, IHLF. (FUNC( I) , I=1,3),DERF( 3) , ( I5QL< I) • 
FORMAT  («  */Tb,*LT  A-*  *FB*5.T22*  *  N—  *  ,  I3,T32,*  I HLF- *  *  I  3 , T40 
6,16) 

IF  ( ISOL(2) .NE.O)  WRITE  (6,20) 

FORMAT  (*  «///) 


1=1,2) 


. 


* 


FI  7 

c 

JO  CONTINUE 
C 

HE TURN 

END 


5 


SUBROUTINE  ST  AT  F2  (  ST  OR  F  A  ) 


FI  8 


C 


c 

c 

c 

c 

c 

c 

c 

c 

£ 

C 

C 

10 

£ 

C 

c 


c 

20 

C 

30 

C 

C 

c 


c 

c 

c 

4  0 

c 


IMPLICIT  REAL  *8  (A-H,0-Z) 

DIMENSION  STQREA( 26. 7 1 . 3 1,  IEMPIJ3  5X».EiiXUt  L*-VIH-6  1 

D I  ME NS I  ON  FV (  7 1  )  •  F X ( 6 ) , AR GY ( 6 ) , ARGX ( 6 ) , V AL ( 6 ) 

EPS=S0-6 

DO  60  K=1 , 3 


D=  1  • 0  1  DO** ( DFLUAT ( K ) /2D0 ) 

DO  60  11=1,3 

X=DSQRT {  IDO/ 1 .0  IDO ) *OFLQAT{  II  ) /I ODO 
DO  60  J  J  =  1 , 3  8 

Y=DSQRT( 1,0100) *DFLOAT ( JJ-1 )/l ODO 
DO  30  1=1,5 


DO  10  J= 1*71 

FY ( J) =STOREA(  I  , J  ,  K ) 

CONTINUE 


CALL  DATSE ( Y , ODO , 0. IDO, FY, 71,1, ARGY , VAL , 6) 
CALL  DAI.  I  (Y  ,  ARGY  ,  VAL,  P,  6,  EPS,  IER  ) 

IJ=1+1 
HD ( I J)=P 


IF  (IER.NE.O)  IER=IER+10 
IF  (IER.NE.O)  WRITE  (6*20)  IER 
IF  ( IER.NE.O)  STOP 

FORMAT  <*  * /T2 , •  I ER= • , I  2/ ) 

CONTINUE 

HD ( 1 )=0D0 

DO  50  J  = l , 6 

J  ARG= 1 ODO  * ( ARGY ( J ) 41 D- 10 ) +1 
DO  40  1=2,6 

FX (  I  )=STGREA(  1-1  , JARG ,K ) 


CONTINUE 


FX ( 1 )=0D0 


J 

'  n 


■ 


FI  9 
c 

CALL  OATSE(  X,  ODO  ,  0 . 1  DO  ,FX  ,6 , 1  •  ARGX  ,  VAL  ,6  ) 


CALL  DALI (X,ARGX,VAL,P,6,EPS,IER) 

£ 


VD ( J) =P 

IF  (IER.NE.O)  IER=IFR+20 

IF  (IER.NE.O)  WRITE  (6,00)  IE  R 

c 

IF  (IER.NE.O)  STOP 

50 

CO NT JNUE 

c 

CALL  DA  TSE ( X  *  OD  0»0,1D0*HQ»6* 1  , ARGX,VAL,6 ) 

CALL  DAL  I ( X , ARGX , VAL, P 1 ,6, EPS, l ER ) 

c 

IF  (IER.NE.O)  ItR=IER+30 

IF  (IER.NE.O)  WRITE  (t>,20)  IE  R 

IF  (IER.NE.O)  STOP 

c 

CALL  DALI ( Y. AR GY , VD , P2 , 6 * EP S, I ER ) 

c 

IF  (IER.NE.O)  IER-IER+40 

IF  (IER.NE.O)  WRITE  (6,20)  l ER 

IF  (IER.NE.O)  STOP 

c 

TEMPI  I I  , JJ,K)=( P1+P2)/ (2D0*D) 

c 

60 

CONTINUE 

c 

DO  70  K=1 ,3 

DO  70  1=1,26 

DO  70  J— 1,71 

c 

STORE A ( I . J ,K )=0D0 

c 

70 

CONTINUE 

c 

DO  80  K=l,3 

DO  80  1=1,3 

DO  80  J=1 ,38 

c 

STOREA ( I , J,K)=TEMP( I , J,K) 

c 

80 

CONTINUE 

c 

RETURN 

END 


c 

c 


MAIN  PROGRAM 


F  20 


c 


c 

c 


c 


c 


c 

c 


c 


c 


c 


c 

c 


c 

c 


IMPLICIT  REAL*8  (A-H,0-Z) 

HIHENSION  Ffj.71.3)  ,FN(3,71.3)  ,FNN(71) _ 

DIMENSION  S  (3,71,3)  ,SN  (3f  71 ,3)  ,SNN(71) 

DIMENSION  PRMT  (5)  ,  FIT  NO  (2)  ,  DERS  (2)  ,AUX  (16,2) 

DIMENSION  STORE  (26,71,3,3)  ,  STORAF  (26,7  1,3)  ,  STORBE  (3  , 7  1 , 3) 
DIMENSION  STORS  (3,71,2,2)  , STOHAS  (26 ,7  f,  2)  , STORES  (  3 , 7 1 , 2 ) 
DIMENSION  SNS  (  3)  , I SOL  ( 2)  , SSP  (71 ,2) 


COMMON  S,SN,SNN,F,FN, FNN , ISOL , NE/FCTD AT/EQN , X,DX ,T, DT,G 
COMMON  / 1 N T D AT/I X , M P S , KT , NS 

EXTERNAL  rCT,OUTP 


10  =  0 


10=  1 

KT  =  1 

KT'=2 

KT=  3 

K  T  =  0 

F.  2 

The  Energy  Equation 

Program 

NS  =  2 

NS=  3 

•  •—  -  •  •  - -----  -  - - - -  •  ■  •-  -  . - 

UE3= 1 . 0 1  DO 

IFX=26 
I F  X  =  3 

IPSR=0 
NDIM=  2 


PPMT  (2)  =  7  DO 
FRMT  (3) =0.  IDO 
PRMT (4) =  5D-  2 
E  R RO  R=S  D-6 


LF  (KT.EC.O)  ITERS=50 

TTERS=  30 
1 1  E  7  S  =  S  0 

DO  30  L  =  1 ,3 

DO  30  J  =  1 . 7  1 

DO  30  1=1,26 

SI 0 RAF  (I  , «J,L)  =0D0 
STORF  (T, J ,  1  , L)  =  0D0 
SI07F  (I .0 .2 , l) =0 DO 

STORF  (1,0, 3, L) =ODO 
IF  (I.GT.3)  GO  TO  10 
ST073F (I, J ,L) =0D0 


■ 


c 

F21 

p 

1  0 

continue 

IF  ( L. SO.  3 )  GO  TO  20 

c 

STORAS  (I, J ,1) =0 DO 

IF  (I.GT.3)  GO  TO  20 

STORS  (I ,  J  ,  1  ,  L)  =01)0 

STORS  (I,J,2,  L)  =  0 D 0 

STORBS  (I.  J.L) =0D0 

c 

20 

CONTINUE 

c 

3  0 

CONTINUE 

c 

DO  40  J= 1 ,71 

c 

FNN  (J) =0D0 

3NN  (J)  =  0D0 

c 

DC  40  1=1,3 

c 

DC  4 n  K= 1 , 3 

c 

f (i, j,k) =  odo 

FN (I, J,K) =0D0 
S (I, J,X) =  0D0 
5  N  (I,,!,  K)  =  0  DO 


c 

40 

CONTINUE 

c 

IF  (KT.  EO.  0)  GO  TO  60 

c 

READ  (7)  STORAS 

c 

DO  50  L=  1 , 2 

DO  50  1=1,3 

DO  50  J=  1 , 7  1 

c 

STORS  (I. J, 1 , L) = STORAS  (I, J,L) 

c 

50 

CONTINUE 

c 

60 

CONTINUE 

c 

IF  fKT. LT. 2. OP. KT. GT. 3)  GO  TO  80 

c 

IF  (KT.BO-2)  READ  (8)  STORES 

IF  (KT. EQ, 3)  READ  (0)  STORES 

c 

DC  7  0  I.=  1,2 

DO  70  1=1.3 

DO  70  J=  1 ,71 

c 

STORS  (I, J,2,L) = STORAS (I,J,L) 

••• 


STORS  (I,  .J,  1rL)  = STORES  (I,  J,L) 
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c 

7  0 

c 

CONTINUE 

8C 

C 

CONTINUE 

IF  (KT. ECU  0 . AND. NS. EQ. 2)  READ  (1)  STORAF 

IF  (KT.EQ. 0. AND. NS. EQ. 3)  READ  (5)  STORAF 

IF  (KT.F0.1)  READ  (2)  STORBF 

IF  (KT.EQ.  2)  READ  H)  STORBF 

c 

IF  (KT. EQ. 3)  READ  (4)  STORBF 

IF  (KT.GE.2)  K=  3 

IF  (KT. EQ.  1)  K  =  2 

IF  (KT. EQ. 0)  K= 1 

c 

c 

DO  100  L=  1 , 3 

DO  100  1=1, IEX 

DO  100  J  =  1 ,71 

STORF (I,»J,K,L)=STORAF  (I, J,L) 

If  rKT.EO.Q1  GO  TO  90 

c 

90 

c 

100 

c 

STORF  (I , J , 1 , L) =STORBF  (I, J, L) 

CONTINUE 

CONTINUE 

c 

IF  (KT. EQ. 2)  READ  (2)  STORBF 

IF  (KT.EQ.3)  READ  (3)  STORBF 

IF  (KT.LT.2)  GO  TO  120 

c 

DO  110  L=  1 . 3 

c 

DC  110  1=1, IEX 

DO  110  J=1 ,71 

STORE  (I, J, 2, L) =STORBF  (I,J,L) 

c 

1  10 

CONTINUE 

c 

1  20 
c 

CONTINUE 

DO  130  1=1,3 

c 

SNS  m  =  0 DO 

c: 

1  30 

c 

CONTINUE 

G=SD0/3D0 

G  =  1 . 4  D  0 

TMACCO= IDO 

c 

E  N  ACCO= 1  DO 

D  X  =  0 .  IDO 

' 


■ 


F23 

DN=0. IDO 

0*1  =  0.0500 

c 

S  N  W  =  0  . 2  DO 

c 

DFNW=0. IDO 

c 

PI=  3.  141 592653589793D0 

DFX=DX 

DTEE=DT 

X  =  0P0 

T  =  0  DO 

C  1=  (2D0-TM ACCO)  /TMACCO 

02  =  200*0/  (G+1D0) * ( 2 DO- EN ACCO) /ENACCO 
'[  F  =  0  D  0 
'I P  =  1  D  0 

TP  =  0. 5 DO 

NF  =  7  1 

C 

DT  =  8  DO*  DTF  E 

IF  (KT.GF.  12)  DT  =  0 . 0  1  DO 

IF  (KT.LF.9)  DT=4D0*DTEE 

Ir  (KT.LF.5)  DT=2D0*DTFE 

IF  (KT.LE.2)  DT= DTEE 

C 

IF  (KT.GE.12)  T=0. 5DO+DFLOAT  (KT-1 1) *DT 

IF  (KT.LE.11)  T  =  0. 3DO+DFLOAT  (KT-9) *DT 

IF  (KT.LF.9)  1  =  0. 1DO+DFLOAT  (KT-5) *DT 

IF  (KT.LE.5)  T=0. 1 DO+DFLOAT (KT-2) *DT 

IF  (KT.LE.2)  T=DFLOAT  (KT)  *DT 

C 

DO  460  I X= 1 , IEX 

c 

POMT  fl)  =0C0 

IF  (KT. EC. 0)  GO  TO  150 

c 

DC  140  K  =  2 , 3 

c 

K  KK  =  K— 1 

c 

DO  140  0=1  ,71 

c 

3 (1, J,K) = STOPS (IX, J, KKK, 1) 

SN  (1, 0,K)  = STOPS  (IX, 0,KKK, 2) 

c 

» 

140 

CONTINUE 

r* 

150 

CONTINUE 

C 

DO  16  0  0=1,71 

c 

FNN (J) =ST0R” (IX , 0 , 1 , 3) 

c 

F (1,0, 1) =  STORF  (IX, J,  1,  1) 

F  (1,Jf  2) =STCPF  (IX, 0,2, 1) 

5 


c 
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F ( 1, J, 3) =5TORF  (IX  ,  J, 3, 1) 

FN  (1, J,  1) =  STORF  (IX,J,  1 ,2) 

FN (1 , J.2) =STORF  (IX.J.2.2) 

c 

FN (1, J,3) =  STO  RF  (IX,  J, 3,2) 

IF  (IX. EQ. 1)  GO  TO  150 

c 

F ( 2,  J  , 1) = STORF (IX- 1 ,J, 1, 1) 

FN(2.J.1)=STORF(TX-1,J.1,2) 

c 

IF  (IX. EQ. 2. OR. IX. EQ. 7. OR. IX. EQ. 14. OR. IX. EQ. 20)  GO  TO  160 

c 

F  (3,  J,  1)  =  STORF  (IX— 2  ,.3,1,1) 

FN (3, J, 1) = STORF  (IX- 2, J, 1 ,2) 

c 

160 

C 

CONTINUE 

C 

IF  (IX. EQ.  1)  GO  TO  180 

IJ=  1 

c 

c 

IF  (IX. SO. 2. OR. IX. EQ. 7. OR. IX. EQ. 14. OR. IX. EQ. 20)  TJ=2 

DO  170  1 1  =  I J  ,  2 

c 

1=4-11 

IJI=I-1 

c 

DO  170  J= 1  ,71 

c 

5  (I, .7,  1 )  =S  (I JI,  J,  1) 

SN  (I,J ,  1)  =  S N  (IJI, J,  1) 

c 

170 

C 

180 

C 

CONTINUE 

CONTINUE 

MDPTCR=0 

LF  =  0 

c 

L  F  1  =  0 

DX=8D0*DEX 

IF  (IX. LE. 19)  DX=4  DO*  DEX 

IF  (IX. LE. 13)  DX=2D0*DEX 

IF  (IX.LF.6)  D  X  =  DEX 

c 

X=X+DX 

ECN=2D0 

IF  (IX. GT.  1 . AND. KT. FO. 0)  E  Q  N  =  3  D  0 

IF  (IX. EQ. 1 . AND. KT. EQ. 1 . OR. KT. EQ. 12)  EQN=22D0 

IF  (IX.GT. 1.AND.KT.EO. I.OR.KT.EO. 12)  EQN=32D0 

I F  (IX. EQ. 1 . AND. KT. GT. 1)  EQN=23D0 

IF  (IX. GT. 1 . AND. KT. GT. 1)  EQN=33DO 

c 


. 


190  CONTINUE 
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C 

I E  S  R  -  0 

_ isol  m  -0 

c 


MPS=1 

IF  (NDPTCR.KQ. 1)  MPS=11 
C 

TF  (KT. EQ.O. AND. IX. GT. 2)  SN W  =  2D0* SN  ( 2 r 1  ,  1 )  - S N  (  3 ,  1  ,  1 ) 

_ IF  (KT.  EG.  11  S N W - S N  (1,1,?) _ 

TF  (KT.GT.1)  SNW=2D0*SN  (1,  1,2)  -SN  (1,  1,3) 

IF  (MPS.EQ.11)  SNW  =  SN  ( 1  ,  1  1 ,  1 ) 

C 

IF  (KT.  EQ. 0. AND. IX. GT.  2)  DSNW=DABS(SN (2 , 1 ,  1)  -SN  (3, 1 ,  1) ) /10D0 
IF  (KT. EQ. 0. AND. IX. GT. 2)  DSNW=DABS(SN (2,1 , 1)  -SN (3,1 , 1)  ) 
_ IF  (KT.  EQ. 0. AND. IX . GT. 2)  PSNW  =  0.01DQ _ 

IP  (KT.F0.1)  DSNW  =  0. 01D0*SN  ( 1, 1 , 2) 

IF  (KT.  GT.  1)  DSNW=DABS  (SN  (1  ,  1 ,2)  -  SN  (1 ,  1  ,  3)  )  /  1 0D9 
IF  (NFS. EQ. 11)  DSNW=0.01DO*SNW 
C 

IF  (MPS.E0.11)  5  N  A  =  S  N  W 


C 


200 

CONTINUE 

c 

IF  (MDPTCR. EQ. 0) 

LF=LF+ 1 

IF  (NDPTCR . EQ.  1) 
ISOL  (2)  =0 

LF 1 =LF 1 + 1 

c 

DIRS  (1)  =0.  3 DO 
DERS  (2)  =0. 7D0 

c 

TF  (MDPTCR.EQ.  1) 

GO  TO  210 

c 

FUNC  (2)  =  S  N  W 

FUNC  (1)  =  TP+('2*FUNC  (2)  /X  +  (C  1  *  *  2-2  D0*C  1 *C  2)  *  (FNN  ( 1 )  /X)  **2-100 


C 

GC  TO  220 
C 

210  CONTINUE 


C _ 

PFNT  (1)  = 1  CO 
F  UNC (2) =  SN  W 
FUNC  (1)  =  S  ( 1 ,  1 1  ,  1) 
C 

220  CONTINUE 


C _ 

SNS  (2) =  FUNC  (2) 

C 

CALL  DHPCG  (P  R  NT , FU  NC, DERS , NDI N , I HLF, FCT  r  OUTP , AU  X) 
C 

IF  (IHLF.GE.11)  WRITE  (6,230)  IHLF 
230  FORMAT  ('  '/TIP, 1 I H  L  F  -  * , 12/) _ 

IF  (IHLF.GE.11)  STOP 

r» 

IF  (DA3S  (FUNC(I)  )  . LT. ERROR)  GO  TO  290 


••• 
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IF  (LF1.15Q.50)  00  TO  290 

TF  (LF.  EQ.  1  .  OR.  ( I. F  1 .  FJ Q .  1  .  AND. M PS.  EQ.  11))  ISOI.  (1)  =  ISOL  (2) 

I F  (IBSR.EQ.2)  GO  TO  270 

IF  KLF.  NR.  1  .  OR.  f LF  1  .  N E.  1  .  ANP.MPS.EO.  1  1 )  )  .  AND.  (ISOL  M)  .  NE. 

c 

1ISOL  (2)  .OR. IBSP.GT.O)  )  GO  TO  240 

SNS  (I SOL  (2) ) =  SNS  (2) 

IF  (ISOL  (2)  .  EQ.  1)  5NW=5NW+DSNW 

IF  (ISOL  (2)  .  HO.  3)  SNW=SN  J-DSN'M 

c 

c 

2  40 
c 

GO  TO  200 

CONTINUE 

I P  SR=  1 

c 

c 

SNS  (ISOL  (2) ) =  SNS  (2) 

II  (DABS  (FUNC  (1) )  . LT. 5D-3)  GO  TO  250 

SN  W  =  (SNS  (  1)  +SNS  (3)  )  /2D0 

IF  (LF+LF 1 . DO. ITERS)  MDPTCB=1 

IF  (LF+LF 1 . F 0. ITERS)  GO  TO  190 

c 

250 

c 

GO  TO  200 

CONTINUE 

I BSP=2 

SF  =  FHNC  ( 1 ) 

c 

SNP=SNW 

S  N  W  =  S  N  S  (  1) 

IF  (ISOL  (2)  .  E0.  1)  SNW=SNS  (3) 

c 

DC  260  J - 1,71 

c 

c 

2  60 
C 

SSP  (J  ,  1 )  =  S  ( 1 ,  J  ,  1 ) 

SSP  (0,2)  =SN  (1,0,1) 

CO  NT  I  N U  F 

c 

270 

C 

GO  TO  200 

CON TIN UF 

SA  =  -FUNC ( 1) /  (SP-FUNC ( 1)  ) 

SP=  1P0-SA 

c 

c 

DC  280  J=1,71 

S  (1 ,0, 1) =SA*SSP  (J, 1)  +  SB*S  (1,0,1) 

SN  (  1 , 0 , 1 )  =  SA  *  SS  P  (0 , 2 )  +  S B *  S  N  (1  , 0 , 1 ) 

c 

280 

c 

790 

CO NTI NUE 

CONTINUE 

' 


\ 
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c 

c 

If  (ME. EQ. 71)  GO  TO  310 

J  S  =  N  F  +•  1 

c 

DO  300  S  ,  7 1 

DFLMDO-P  FLOAT  (71-J)  /DFLOAT  (71~NE}_ 

3 (1,  J,1)=DEI*S  (1,J,1) 

SN  ( 1,  J, 1 ) =DEL*SN  ( 1 , J,  1) 

SNN  (J)  =DEL*SNN  (J) 

c 

3  OC 
C 

3  1  0 
C 

CONTINUE 

CONTINUE 

IF  (MFS.EQ.I)  GO  TO  330 

c 

c 

DFLSN= (SN (1,11,1) -SNA)  /2DQ 

DC  3  20  J  =  2 , 2  0 

c 

El  A  FLOAT  (J-1)  /10DQ 

c 

3  20 

IF  (J.LE.11)  SN  (1,J, 1)  =  SN  (1 ,0, 1) +ETA*DELSN 

IF  (I.OT.  1  1)  S  N ( 1 , J, 1) =S  N (1 ,  J , 1 ) - (2 DO- FT A) * PELSN 

CUN  ]? I N u i: 

U 

330 

CONTINUE 

C 

DO  3  40  J  = 1  , K E 

c 

s (1 ,1, 1 ) =S  (1 ,  J, 1)  +1D0 

c 

340 

CONTINUE 

C 

FSETA=DN*DFI,OAT  (NE-1) 

C 

3  50 
C 

If  (IX. EQ.  1)  WRITE  (6,35)) 

FORMAT  (M‘) 

36  0 

I  S=T 

If  (KT. EC. 0. AND. NS. EQ. 3)  T=10D0 

WRITE  (6,360)  X, FSETA , T 

T-TS 

FORMAT  ( ' 6 •  ,T  10, *  X= '  ,F4. 1 ,T20 , ' ETA  (FREE  STREAM)  =  ’ , F3.  1 , 

1T45. *T= 1 ,F6. 4) 

C 

370 

WRITE  (6,370) 

FORMAT  (•  '///T 10, ’UNSTEADY  STATE  (WITH  SLIP)  ENTHALPY  DISTRIBUTED 

IN:  V) 

C 

WRITE  (6.380)  ( S  (  1 , J ,  1 )  ,  U  -  1  ,  N  E ) 

3  80 

C 

FORMAT  ('  1 ,T2, 5F23. 5) 

WRITE  (6,300) 

■.  •  *5 
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3  90 

FORMAT  ('  ’ ///T 10  ,  *  UNSTEADY  STATE  (WITH  SLIP)  S-PRIME  DISTRIBUTION 

c* 

1  :  '  /) 

WRITE  (6.3301  IS N  l  1  ,  J  .  1 1  .  1  .  NE1 

c 

c 

LFS=LF+LF 1 

WRITE  ((>,400)  LP,LF1,LPS 

4  0  0 

FORMAT  ('  ' /////? 10, ’CONVERGED  IN  ' , IP , 1 + '  ,  12 , '  =  '  f T2 ,  '  ITE 

1  1  ) 

RATIONS. 

c 

WRITE  (6,350) 

c 

DO  4  10  J  =  1  ,  N  H 

c 

S(1.J.1)=S(1iJ,11-  IDO 

c 

4  10 

CONTINUE 

c 

IF  (KT.EO.O)  GO  TO  430 

c 

JT=2 

IF  (KT. EO. 0. OR. KT. EQ. 2. OR. KT. EQ. 5. OR. KT. EQ . 9. OR . KT . EQ . 1 1 ) 

JI=1 

c 

n 

DO  420  K=1,JI 

DC  4  20  J -  1 , 7  1 

c 

STORS  (TX,J,K,1)=S  ( 1  ,  J  ,  K ) 

STORS  (IX  ,J,K  ,2)  =SN  (1  ,  J,K) 

c 

4  20 

CONTINUE 

C 

4  30 

CONTINUE 

C 

IF  (KT.GT.O)  GO  TO  450 

C 

DO  440  J  =  1 , 7  1 

c 

ST  ORA  S  (IX  .J  .  1)  =S  (  1 ,  J,  1) 

STORAS (IX, J, 2) =  5N  (  1  , J,  1) 

c 

440 

C 

CONTINUE 

4  50 

CONTINUE 

C 

46  0 

CONTINUE 

C 

IF  (IC.EQ.O)  GO  TO  490 

IF  (KT.NF.O)  GO  TO  4  70 

c 

IF  (1MS.E0.2)  WRITE  (7)  STORAS 

IF  (NS. FQ.  3)  WRITE  (11)  STORAS 

c 

GO  TO  490 

' 
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C 

4  70  CO  NTT  NTJ* 
C 


no  goo  r=i.3 


no 

4  80  L=  1  ,2 

DO 

480  J=1,71 

c 

STORRS  (I  ,v7 ,1)  =  STOPS  (I  ,  J,  1  ,  L) 

c 

480 

CONTINUE 

C 

IF 

(KT.  EQ.  1) 

WRITE 

(8) 

S  rOHBS 

IF 

(KT.  EQ.  2) 

WRITE 

(8) 

S  TORES 

IF 

(KT.  EQ.  I) 

WRITE 

(10) 

STORBS 

C 

490 

CONTINUE 

c 


S  T  O  P 
END 


c 


c 


c 

c 


c 


c 

c 

c 


SUBROUTINE  FCT  (ETA  ,  FTJNC,  DERS) 

IMPLICIT  REALMS  (A-H,0-Z) 

iIlMENBION  F  (3.71.3)  ,FK(3,71,3)  ,FNN(71) 

DIMENSION  S  (3,71 , 3)  ,SN  (3,71 , 3)  ,  SNN  (71) 
DIMENSION  FT  (61 ,2)  ,  ARC  (20)  ,  VAL  (40) 
DIMENSION  FUNC  (2)  ,  DERS(2) 

CO U MON  S  ,  S  N  ,  S N N  ,  F  ,  FN  ,  I*  N N 

COMMON  /FCTD AT/EQN , X , D  X  ,  T  ,  D  T  ,  G _ 

COMMON  /I NT DAT/ IX , MPS ,KT, NS 

ERR=5D-6 


N  1  =  M  PS  +  60 
N  2  -  20 

N  3  =  6  1 
EFS=ERR 


B  G  N  =  0  D  0 

IP  (MPS.  E  0 •  1 1 )  BGN= 1  DO 


DFRS  (1 ) =FUNC  (2) 

I N  TE  3  ?= 1 

DC  10  J=  M  P S , 7  1 
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N=  J 

c 

IT ACHK=0 .  1  DO  *D FLOAT  (N-1) 

D  F  T  A  =  ET  A-  E T  A  C  H  K 
C 

_ IF  (DABS  (BETA)  . LE. 5P- 1 0)  IMTERP=0 

IF  (INTERP. EQ. 0)  ET  A=  ETACFK 
IF  (INTERP .  F'Q.  0)  GO  TO  20 
IF  (DETA. LE. ODO)  GO  TO  20 
C 

10  CONTINUE 

C _ 

20  CONTINUE 
C 

FAXT=F(1,N,1) 

FNAXT=FN  (  1  ,N , 1) 

C 

_ FM 1X  =  F  (2.N  .  11 _ 

F  M  2  X=  F  ( 3  ,  N  ,  1 ) 

C 

FM1T=F ( 1 ,N ,2) 

FM2T=F ( 1 , N , 3 ) 

C 

_ SM  1X=S  (2.N  .  1) _ 

SM2X=S(  3,1, 1) 

C 

S  M 1 T=  S  (  1,1,2) 


I  5 


. 


. 


. 


' 
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Sr2T-S ( 1 ,  N  ,  3 ) 

c 

II  (INTIS! T  .  ?v.  0)  GO  TO  110 

C _ 

1K  =  0 

DC  30  J  =  M  F  S  ,  N  1 
C 

J  J  =  J-!TPS  + 1 

FT  (JJ  ,  1)  =F  (  1,J,  1) 

FT  (JJ.2I=PN(1,J,1) 

C 

30  CONTINUE 
C 

I K  =  T  K  +  1 

CALL  DATSF  (  ETA  ,BGN  ,  0.  IDO  ,  FT,  N3 , 2  ,  AEG  ,  V  AL,  N2) 

_ CALL  PAH I  ( FT  A , ARG , V  A  L , P , N2,BFS,IEP) _ 

F  A  XT=P 
C 

IF  (IFP.GT.1)  WRTTE  (6,40)  IER 

40  FORMAT  f1  ’/T  1 0 , 1 1 ER= ' , 1 1 , '  FOR  FAXT. '/) 

II  (ILR.F0.1)  EPS=1 0D0*EPS 

_ IF  (IER. EO.  1  . AND. IK. IT. 5)  GO  TO  30 _ 

II  (IFR.IC.1)  WRITE  (6,40)  IFF 

IF  (IFR.OT.O)  STOP 

EPS=  FFR 
C 

I  K  =  0 

_ DO  60  J  =  f1PS.Nl _ 

C 

J J=J-MPS+ 1 

FT  (JJ,  1)=FN(1,JX1) 

FT  (JJ  ,2)  -  FN N  ( J) 

C 

50  CONTINUE _ 

C 

I K  =  I K  +  1 

CALL  DATSF  (  FT A , flGN , 0 ,  1  DO , FT , N3 , 2 , AEG , V  4  L , N  2 ) 

CALL  DA  HI  ( ET  A  ,  ARC  ,  V  A  I  , F,N2 , EES , TER) 

F  N  A  X  T = P 

C _ 

II  (IFP.GT.1)  WRITE  (6,60)  IFF 

60  FORMAT  (’  ' /T  1  0 , ' IER= ' , 1 1 , *  FOR  FNAXT.  V) 

IF  (IER.FQ.1)  FPS=10p0*EPS 
IF  (IER. EQ.  1  . AND. IK. LT. 5)  GO  TO  50 
IF  (IEP.E0.1)  WRITE  (6,6  0)  IER 

_ IF  (IFR.OT.O)  STOP _ 

I PS=FPR 
C 

IK  =  0 

DO  7Q  J  =  N F S , N 1 
C 

_ j  j  -  j  -  n  p  .s  + 1 _ 

FT  (JJ,  1)  r:F  (2,  J,  1) 

FT  ( J  J  ,  2  )  =■  F  N  (2,  J,1) 

C 


585 


<■ .  *tv  i 


. 

. 


3 


7  0 
C 
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CONTINUE 

I K  =  I K  +  1 

CALL  DATSL  (ETA.BGN.O.  1  00  ,  gT  .  N  3 . 2  .  A  RQ  ,  V  A  L  .  N2  ) 

c 

CALL  DA HI  (ETA, ARG,VAL, P, N2, ELS, IER) 

FM  1  V=P 

HO 

II'  (IER. (IT.  1)  WHITE  (6, HO)  IER 

FORMAT  ('  '/TIC, * IER=' , 11 , »  FOR  FM1X-'/) 

IF  (IER.E0.1)  EPS=10D0*EPS 

C 

TF  (IER. EQ.  1.ANP.TK.LT.5)  00  TO  70 

IF  (TEH. EC. 1)  WRITE  (6,80)  IER 

IF  (IER.GT.O)  STOP 

EPS=FRR 

I  K  =  0 

C 

DO  90  J  =  M P S  ,  N 1 

J J= J-M PS+  1 

FT  (JJ, 1) =S  (2,  J,  1) 

FT  (JJ,  2)  =  SH  (2,  J,  1) 

c 

9  0 
C 

CONTINUE 

rK=iK+i 

CALL  DAIS!  ( FT  A , RGN ,  ) .  100, FT, N  3 , 2 , A  R  0 , 7  A  L , N  2 ) 

CALL  UAHI  (ETA ,ARG,VAL, P , N2 , E PS , IER) 

SHU  =  P 

c 

100 

IF  (IER. 01.1)  WRITE  (6,100)  IER 

FORMAT  ('  ’/T10, ' IEH=’  ,11  , *  FOR  SM1X.'/) 

IF  (I ER . EO . 1 )  EPS= 1 OD 0* EPS 

IF  (IER. EQ. 1 . AND. IK. LT. 5)  GO  TO  90 

IF  (IER. EC. 1)  WRITE  (6,100)  TER 

c 

1  10 

c 

IF  (IER.GT.O)  STOP 

F  P  3  =  E  R  R 

CONTINUE 

IF  (EON. NE. 2 DO)  GO  TO  120 

c 

DERS  (2)  =-  (FA  XT+ETA)  *F!JNC  (2)  -X*  (  (FAXT-FM1  X)  *FUNC  (?)  -  (FUNC  (1)  -30  IX)  * 

1  (F  NAXT+ IDO)) /DX 

c 

RETURN 

c 

120 

c 

CONTINUE 

IF  (IKTE3E.EQ.O.OR. EQN. EO.3B0)  GO  TO  1?) 

c 

TK  =  0 

DO  1  30  J  =  M  P  3  ,  N  1 

c 

J  J  =  J  -  M  P  3  +  1 

FI  (JJ,1)=F.(1,Jr?) 

5 


' 


■ 


I 

D 


c 

1  3  0 

F'l  (JJ  ,  2.)  =FN(1fJ,2) 
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CONTINUE 

c 

c 

IK  =  IK  +  1 

CALL  PATSE  (ETA  r  RGN, 0. 1  DO , FT , N3 , 2 , ARG r V AL , N 2 ) 

CALL  DA  HI  (ETA  ,  A  F  G  ,  V  A  L  , P , 3  2, EPS, TER) 

FM  1  T=P 

IF  (IFR.GT. 1)  WRITE  (6.140)  IER 

140 

FORMAT  ('  */Tl0,  ' IER='  ,11 ,  '  FOR  EMIT.'/) 

IF  (IER. ECU  1)  EPS=10D0*EPS 

IF  (TFH. FO. 1 . AND. IK. LT. 5)  GO  TO  130 

IF  (IER. EC. 1)  WRITE  (6,140)  IER 

IF  (IFR.GT. 0)  STOP 

FFS=EP R 

c 

IK=0 

c 

DC  ISO  J=HFS,N1 

JJ=J-KPS+ 1 

FT  (JJ,  1)  -S  (1,  J,2) 

r 

FT  ( JJ , 2 ) =S N  ( 1  , J  ,  2 ) 

ISO 

CONTINUE 

c 

I K  =  I K  +  1 

CALL  PATSF  (ETA,BGN,0.  1  D 0  ,  FT  ,  N3 , 2  ,  A RG  ,  V  A  L  ,  N2 ) 

CALL  DAH I  ( FT A , A  RG , V  A  L , F,  N  2  ,  F  F  5  ,  I  E  R ) 

S»1T=P 


c 

IF 

(IFR.GT.  1)  WPITF  (6,160) 

IF  R 

160 

FORMAT  ('  * /T 10,' IE R=’, 11,' 

FOR 

SHIT.  V) 

IE 

(IER. EO • 1 )  EPS=10D0*EPS 

IF 

(IFR.EQ.  1. AND.IK.LT. 5)  GO 

TO 

150 

IF 

(IER. EQ. 1)  WRITE  (6,160) 

IER 

IF 

(IFR.GT. 0)  STOP 

EPS 

=  E  F  R 

c 

170 

CONTINUE 

c 

IE 

(EON . N  F . 2  2  DO )  GO  TO  130 

c 

1)FRS  (?)  =-  (FAXT+ETA)  *FUNC  (2)  -X*  (  (FAXT-EM  IX)  *FIJNC  (2)  -  (FUNC  (1)  -S M  IX)  * 
1  (FNAXT+1DP)  ) /DX+2D0*  (1D0-T)  **2*  {FUNC  (1)  -SMlT) /DT-2DO*T*  (1DO-T)  *  (  (F 
2 1)  N  C  (1)  -  S  FIT)  *  (FNAXT+ IDO)  -  (FAXT-FMlT)  *FUMC  (2)  )  /DT 


C 


RETURN 

c 

1  80 

CONTINUE 

C 

IF  ( I N T F, F; T  .  1- 

0.0)  GO  TO  230 

c 

IK=0 

■  • 


' 


c 
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JJ=J-MPS+  1 

FT  {J.J,  1)  =F  (3,0,  1) 

PI  (JJ.2)  =  FN  ( 3  .  J  .  1) 

c 

100 

c 

CONTINUE 

I  K=  T  K  +  1 

CALL  DATSF,  (ETA  ,BGN,  0.  IDO  r  FT,  N3 ,2  ,  ARG  , VAL, N2) 

CALL  DA  HI  (FT  A.ARG.VAI..P.  N2.EFS.IER) 

c 

FM2y=P 

200 

IP  (IEB.GT.1)  WRITE  (6,200)  IER 

FORMAT  ('  1 / T 1 0 , *  I E  R  =  ' ,11,  '  FOR  FM2X.'/) 

IF  (TER. EC. 1)  EPS=10DO*EPS 

IF  (IEP.FO.  1. AND.IK.LT. 5)  GO  TO  190 

c 

TF  (IFF. F  0 . 1 )  WRITE  (6,210)  IER 

IF  ( I f; P  .  G 1 . 0 )  STOP 

F.PS  =  F,RR 

I K  =  0 

DO  210  J=MPS.N1 

c 

J J  =  J-MPS  +  1 

FT  (JJ,  1)  =S  (3,0,1)  .  _ 

FI  (JJ, 2) =SN (3, J,  1) 

c 

210 

CONTINUE 

C 

I  K  =  IK  +  1 

C A I L  DATSE  ( P T A , BGN , 0 .  IDO, FT , N3 , 2 , APG , V A L , N 2) 

CALL  D  A  B I  ( E T  A  ,  A R G  ,  V A  I. ,  P ,  N2  ,  pfp S  ,  IEP )  ‘ 

SM  2X=  P 

C 

220 

IF  (IFR.GT.1)  WRITE  (6,220)  IER 

FORMAT  ('  '  /T  1  0  ,  '  I£.R  =  1  ,  T  1  ,  '  FOR  SH2X.  »/) 

IF  (IER.EQ.1)  EPS=10D0*EPS 

IF  (IER. EC. 1.AND.JK.LT. 5)  GO  TO  210 

IF  (IER.EQ.1)  WRITE  (6,220)  IER 

IF  (IFR.GT.O)  STOP 

C 

230 

C 

EFS-EPR 

CONTINUE  _ _ _ _ 

IP  (EQN.NF.3D0)  GO  TO  240 

C 

C 

DFR5  (2)  =-  (FAXT  +  ETA)  *FUNC  (2)  -X*  (  (1 . 6  DO  *  F  A  XT- 2D0  *F  M  1  X  +  0 . 5  DO  *  F  M  2X)  *ftj 
INC  (2)  -  (1. 5D0*FUNC  (1)  -2D0*SM1X+0.5P0*SM2X)  *(FNAXT+1D0)  )  /DX 

P FT URN 

c 

24  0 

CONTINUE 

c 


c 


IP  (ECN.NP.32E0)  GO  TO  230 


. 
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Di'RS  (2)  =-  (EAXT+ETA)  *FIJNC  (7)  -X*  (  (1 . 5D0*  FA  XT- 2D0*FM  1  X  +  0 .5D0  +  FH2X)  *FU 
INC  (2)  -  (  1.  5  DO*  FIT  NC  ( 1)  -2  DO*  SMI  X  +  0.  5D0*SM2X)  *  (FNAXT  +  1D0)  )  /DX  +  2D0*  { IDO 

2- -i)  **2*  (FUNC  (1)  -  SM1T)  /0T-2D0*T*  (1D0-T)  *  (  (FUNC  (1)  -SNIT)  *  (FNAXT+  IDO) 

3-  (F AXI-FM  1 1)  *  F U  NC  (2\  )  /DT 

c 

BFTURN 

c 

2  50 

c: 

CONTINUE 

f 

IF  fINTEEP. EO. 0)  GO  TO  300 

c 

c 

I K  =  0 

DO  260  J=MPS,N1 

DJ  =  J-MPS-*-1 

FI  (JJ,1)  =F  (1fJ(3) 

c 

25  0 
C 

FT  (JJ,2)=FN  ( 1  ,  J  ,  3) 

CONTINUE 

I K  =  I K  +  1 

CALL  DAT  Si-,  ( E  T  A  ,  b  G  N  «  0  .  1 D  0  ,  F  T  ,  N  3 , 2  ,  AF.G  ,  V  A  L  ,  N2 ) 

c 

CALL  D  A  H I  (  ET  A  r  AF.G  ,  V  A  L ,  P ,  N2  ,  EPS  ,  I ER) 

F  M  2  T  =  P 

270 

IF  ( I F  K . G  T .  1 )  WRITE  (6,270)  IFP 

FORMAT  (•  »/T  10, ' IEP= '  , 11 ,  ’  FOR  FM2T. '/) 

IF  (IFR.E0.1)  EPS=10D0*E?S 

C 

IF  (IFF. EO. 1. AND.IK.LT. 5)  GO  TO  260 

IF  (IEh. E0.1)  WRITE  (6,270)  IE  R 

IF  (ISR.GT.O)  STOP 

EPS=EJRH 

I K  =  0 

c 

DO  280  J=MPS,N1 

JJ=J-MPS+  1 

FT  (JJ,  1) -S  0,J,3) 

FT  (JJ,  2)  =  S  N  (  1 ,  J,  3) 

c 

280 

c 

CONTINUE 

i  ;<  - 1  k  + 1 

CALL  3 A T SE  ( ^n'  A  ,  hGN ,  0  .  1  DO  ,  FT  ,  N  3 , 2  ,  APG  ,  V  A L  ,  N  2 ) 

C  A I.  L  DA  H  I  {  F‘r'  A  ,  A  K  G  ,  VA  L  ,  1  ,  N  2  ,  E  P  S  ,  I E  R) 

3  w  2  T  =  P 

c 

290 

IF  (IFR.GT.1)  WRITE  (6,2)0)  IEP 

FORMAT  ('  '/T 1 0 , 1 TER= ’ , I 1 , '  FOR  SM2T. '/) 

IF  (IER.E0.1)  EP  S= 1 0D0*E  PS 

TF  ( I E  R . E 0 •  1. AND.IK.LT. 5)  GO  TO  280 

IF  (IER.EP.1)  WRITE  (6,2)0)  TER 

TF  (TER.GT.O)  STOP 

E  P  S  =  E  R  R 

c 


■ 


\ 


\ 


30  0  CONTINUE 
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C 

IF  (EON . N E. 2  IDO)  GO  TO  310 

£ _ 

DERS  (2)  =-  (FA  XT* ETA)  *FUNC  (?)  -X*  (  (FAXT-FM  IX)  *FUNC  (2)  -  (FUNC  (  1)  -SU  IX)  * 

1  (FNAXT  +  IDO)  )  /  DX  +  2D0*  ( 1D0-T)  **2*  (1  .'r>n0*FUNC  (1 )  -  2D0*3M  1T  +  0. 5nO*5U2f) 

2/ET-2D0*T*  (1D0-T)  *  (  (1. 5D0*FUNC  (1)  -2D0*SU  1T+  0 . 5D0*  S  * 2T)  *  (FNAXT+1D0) 
3-  ( 1 .  SD0*FAXT-2D0*Ff11T+0.  ‘3DQ  +  FM2T)  *FUNC  (?)  )  /DT 
C 

_ RETURN _ 

C 

310  CONTINUE 
C 

DERS  (?)  =-  (FAXT  +  ETA) *FUNC  (2)  -X* (  (1 . 5  DO  *  F  A  XT- 2D0  *F  M 1 X  +  0 . 5D0*FM2X)  *FU 
INC  (2)  -  (  1. 5  DO  *  FUNC  (1)  -2D0*S1''!1X  +  0.5D0*SI'12X)  *  (INAXT  +  1D0)  )  /PX  +  2D0*  (100 
_ 2-T) **?*(1,5P0*FUNC(1) -2D0*SM1T+Q. 5PQ*SM2T) /DT~2D0*T*  (1P0-T)  *  (  ( 1 . 5D 

3  0*  FUNC  (1)-2P0*SM1T+0.5D0*SF2T)*  (FNAXT+  100)  -  ( 1. 5D0*FAXT-?D0*FN1,J>  0. 

4  5D0*FM2T) *FUNC (2)  ) /DT 
C 

RETURN 
£  V  D 


S UR ROUT  INF  OUTP  (ETA , FU NC , P P R S , I HLF , N  D I  M! , PPMT) 
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C 

IMPLICIT  IEAL*fl  (A-H,0-Z) 

DIMENSION  F  (3.71 .3)  , F N  (  3 .  71.3)  ,FWN  (7 1) 


L  J  M 

ENSION  S (3,71,3)  ,  S  N  ( 3  , 

71,3)  ,  SNN  (71) 

DIMENSION  FRMT  (5)  , FUNC  (2) 

, DEES  (2)  , ISOL  (2) 

c 

COMMON  S  ,  S  N  , 

SNN,F,FN,  INN, 

IS Cl  ,  N 

COMMON  /INTPAT/IX,MPS,KT, 

NS 

c 

IE 

(KT. NF. 0) 

GO  TO  2 

c 

IF 

(IX.  EQ.  1) 

N  E  =  4  2 

IF 

(IX. EQ. 2) 

N  E  -  4  4 

IF 

(IX. EQ. 3) 

N  E  =  4  5 

IF 

(NS. EG. 2) 

NE=39 

c 

2 

CONTINUE 

' 

c 

IF 

(KT. NF.  If 

GO  TO  4 

c 

IF' 

(IX. EG.  1) 

N  E-4  0 

IF 

(IX.  EQ. 2) 

NE  =  4  1 

IF 

(IX. EQ. 3) 

N  E  =  4  2 

c 

4 

CONTINUE 

c 

IF 

(KT. NE. 2) 

GO  TO  6 

c 

IF 

(IX. EQ.  1) 

NE  =  4  1 

IF 

(IX.  EQ. 2) 

NE  =  4  2 

IF 

(IX.  FQ.  3) 

NE=42 

c 

f, 

CONTINUE 

c 

IF 

(KT. NE. 3) 

GO  TO  7 

c 

TF 

(IX. EQ.  1) 

NE=4  2 

TF 

(IX. EQ. 2) 

NE“4  3 

IF 

(IX.  EQ. 3) 

NE-4  4 

c 

7 

C  C  N 

TIN  HE 

c 

NX  = 

0 

NX  - 

3 

xi  V  - 
»  ■ 

2 

NX  = 

1 

N  X  = 

4 

N  X  = 

5 

c 

INTERP= 1 

c 

DC 

10  J  =  1 , 7  1 

C 


N  =  J 
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C 

ETACHK=0.  1D0+DFLOAT  (N-1) 

DFT  A=  ETA- ETA  CHK 

£ _ 

IF  (DABS  (DETA)  .LE.  5D-1Q)  INTERP=0 
IF  (INTERP.EQ. 0}  ET  A  =  ETACHK 
IF  (INTFPF. EQ. 0)  GO  TO  20 
IF  (DETA . LE. ODO)  GO  TO  20 
C 

10  CONTINUE 

C 

20  CONTINUF 
C 

IC  =  1 

IF  (IX.GT.NX)  .10=1 
10  =  0 

C 

IF  (INTERP . ME. 0)  RETURN 
C 

S  (1,N,  1)=FUKC(1) 

S N  (  1 ,  N  ,  1 )  =  FU  NC  (2) 

SNN  (N)  =  PERS  (21 

C 

IF  (FUNC(1)  .GT.ODO)  ISOL(2)=3 
IF  (FUNC  (1)  .  IE.  ODO)  IS0L(2)=1 
IF  (DABS  (DEPS  (2)  )  .GT.  1D2)  PRMT(5)=1D0 

IF  (W.EQ.NE)  PPMT  (5)  =  1  DO 

C _ 

IF  (10. EO- 0)  GO  TO  50 
C 

WRITE  (6,3  0)  ETA , N , I  HI F,  (FUNC  (T)  ,1=1,2)  , DER  S  (2)  €  (TSOL(I)  ,1=1,2) 

30  FORMAT  (’  * /T5, * ETA= • , F8 . 5 , T22 , ' N= 1 ,I3,T32, *IHLF=» , 13 ,T4 0 , 3 D20 . 0 , T 

16,16) 

C _ 

IF  (N.EQ.RE)  WRITE  (  6  ,  U  0) 

40  FORMAT  ('  '///) 

C  . . ______ . „ . . . 

50  CONTINUE 
C 

EETUPN 


END 


■ 


